Monday, February 17, 2014

Efficient Trading Option Pricing in Scala

Target audience: Intermediate
Estimated reading time: 5'

Ever pondered how the pricing of an option on a future or index works? Curious about the underlying formula?
This article offers a succinct summary of the Black-Scholes formula and discusses how it's implemented in Scala, utilizing memoization to enhance execution speed.


Table of contents
      Formulation
      Memoization
Follow me on LinkedIn

What you will learn: A fast, memoized implementation of option pricing equation of Black-Scholes.

Notes
  • 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.
  • Familiarity with option trading [ref 1] is not required to understand the topic
  • Environments: Scala 2.12.10, JDK 11

Overview

You don't need to be a seasoned professional in option trading or asset valuation to grasp the essence of Black-Scholes. 

The Black-Scholes model is a method employed to calculate the fair or theoretical value of a call or put option. This calculation is based on six key factors: the volatility, the type of option, the price of the underlying stock, the duration until expiration, the strike price, and the risk-free interest rate.

The Black-Scholes-Merton equations, a notable set of stochastic partial differential formulas, are renowned for pricing options based on underlying security's price, volatility, and expiration time [ref 2]. Functional languages like Scala give developers the flexibility to experiment with and test various volatility and time decay models [ref 3].

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 produces the following solution for a call option: \[(1)\ \ C(t) = N(d)S - N(d')Ke^{-r(T-t)}\,\,with\\N(x)=\frac{1}{\sqrt{2 \pi }}\int_{-\infty }^{x}e^{-\frac{t^{2}}{2}}dt\\(2)\ \ d = \frac{1}{\sigma \sqrt{T-t}}[log(\frac{S}{K}) + (r+\frac{\sigma ^{2}}{2})(T-t)]\\(3)\ \ d^{'}=d-\sigma \sqrt{T-t}\] K is the strike price for the underlying security. The price of a put option is computed as \[P(t)=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

Calculating the price of an option requires extensive CPU usage. The approach includes several steps:
  1. Establish and encapsulate the parameters used in the Black-Scholes formula.
  2. Carry out the process to calculate the price of a call option.
  3. Identify the various computational elements involved in the Black-Scholes formula.
  4. Create a system for storing temporary results.
  5. Combine these interim outcomes to finalize the call price.
  6. Use Monte Carlo Simulation for determining the security's price.

Formulation

First let's define a data class, BSParams to encapsulate the parameters of the Black-Scholes formula:
  • S Current price of the underlying security
  • K Strike price
  • r Annual compound risk-free interest rate
  • sigma Volatility of the price of the underlying
  • T Expiration date for the option
case class BSParams(
   S: Double, 
   K: Double, 
   r: Double, 
   sigma: Double, 
   T: Double)


The Scala code snippet below implements the pricing of a call option for a security S and reflect the mathematical equations defined in the previous paragraph.

def callPrice(p: BSParams, t: Double): Double = {
    import Math._

   // Implement the Black-Scholes stochastic equation,
   val sigmaT = p.sigma * sqrt(p.T - t)

   // Formula (2)
   val d = (log(p.S/p.K)+ (p.r + 0.5*sigma*sigma)*(p.T-t)) /sigmaT 

   // Formula (3)
val d' = d1 - sigmaT
   
val gauss = new Gaussian
   // Computation of the call C(t) as per formula. (1)
   p.S*gauss.value(d) - p.K*exp(-p.r*p.T)*gauss.value(d')
}

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. 

Memoization

Memoization, in the realm of computing, is a strategy for optimizing program performance. It involves caching the outcomes of costly function calls to pure functions. This way, when inputs are repeated, the program can quickly retrieve the stored result instead of recalculating it.

The aim is to reduce computational expenses by incorporating the algorithm into a workflow. In this approach, steps within the workflow that remain unchanged by the alteration of a single variable's value are not executed. For example, a change in the underlying security's price, denoted as S, necessitates recalculating only a specific part of the workflow, namely \[S=> log(S/X)\] The Black-Scholes algorithm is divided into five distinct internal computational stages:

object BSParams{
   val fLog = (p: BSParams) => Math.log(p.r / p.T)
   val fMul = (p: BSParams) => p.r * p.T
   val fPoly = (p: BSParams) => 0.5 * p.sigma * p.sigma * p.T
   val fExp = (p: BSParams) => -p.K * Math.exp(-p.r * p.T)
}

The 4 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, K, r, sigma or T is changed.
The computational state is completely described by
  • Parameters of Black-Scholes equations S, K, r, sigma, T
  • Partial or intermediate results {bs1 to bs5}
Next, we need to wrap the computation and update of intermediate results of the Black-Scholes formula.

case class PartialResult(bs1: Double, bs2: Double, bs3: Double, bs4: Double, bs5: Double){
    import BSParams._

    lazy val d1: Double = (bs1 + bs2 + bs3) / bs4

    def f1(p: BSParams): PartialResult = this.copy(bs2 = fLog(p))

    def f2(p: BSParams): PartialResult = this.copy(bs2 = fMul(p), bs5 = fExp(p))

    def f3(p: BSParams): PartialResult = this.copy(bs3 = fPoly(p), bs5 = fExp(p))
 }

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

class State(p: BSParams, pr: PartialResult){

    def setS(newS: Double): State = new State(p.copy(S = newS), pr)

    def setR(newR: Double): State = new State(p.copy(r = newR), pr)

    def setSigma(newSigma: Double): State = 
                  new State(p.copy(sigma = newSigma), pr)

    lazy val call: Double = {
       import org.apache.commons.math3.analysis.function.Gaussian

       val gauss = new Gaussian
       p.S * gauss.value(pr.d1) -pr.bs5 * gauss.value(pr.d1 - pr.bs4)
    }
}

Each of the methods setS, setR, setSigma updates one of the Black-Scholes equation parameters that automatically triggers a selective re-computation of some of the intermediate results.
The price of the call is a lazy value computed only once (immutable state).  The update of the state of the Black-Scholes computation is rather straight forward.
The following code snippet illustrates the immutable update of sigma in the price of a call 

// 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
log.info(callValue)
// Update the value of sigma val state1 = state0.setSigma(0.14) val newCallValue = state1.call
log.info(newCallValue)

Monte-Carlo Simulation

A Monte Carlo simulation is a method employed to estimate the various possible outcomes of a process that is challenging to predict owing to the presence of random elements. This technique is instrumental in assessing the effects of risk and uncertainty.
The Monte Carlo simulation finds application in an array of disciplines, such as finance, business, physics, and engineering, among others. It's often used to solve a variety of problems and is sometimes known as a multiple probability 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.
The price of security S(t) at time t is defined as \[(4)\ \ S_{t}=S_{0}+e^{(r - \frac{\sigma^{2}}{2})t+\mathbb{N}(0,1)\sigma\sqrt{t}}\]
The solution S(t) can be approximated using an iterative or recursive process. \[(5)\ \ S_{t}=S_{t-1} + e^{(r - \frac{\sigma^{2}}{2}).dt+\mathbb{N}(0, 1)\sigma\sqrt{dt}}\] 
This particular implementation uses a fast tail recursion to simulate the computation of the underlying security at any given time t. The sampling times are normalized as integers with dt =1 for a faster simulation.

 import Random._

def monteCarlo(p: BSParams, s0: Double): List[Double] = {
     // Implementation of equation (4)
   val initialValue = s0*exp((p.r-0.5*p.sigma*p.sigma)+p.sigma*nextGaussian)
   monteCarlo(p, 0, List[Double](initialValue))
}

  
@tailrec
private def monteCarlo(p: BSParams, t: Int, values: List[Double]): List[Double] = {
   import Math._
    
   if(t >= p.T)
     values.reverse
   else {
        // Implementation of equation (5)
val st = values.last * exp((p.r-0.5*p.sigma*p.sigma)+p.sigma*nextGaussian) monteCarlo(p, t+1, st :: values) } }

The implementation of the Monte Carlo simulation relies on the fast Scala tail recursion [ref 4].

Evaluation

Let's predict the price of security using the Monte Carlo method for 100 units of time.

val S0 = 80.0
val r = 0.06
val sigma = 0.2
val K = 0.9
val bsParams = BSParams(S0,  r, K, sigma, 100)
val pricePrediction = monteCarlo(bsParams, S0)

The plot below illustrates the impact of the risk free of interest rate, r (6% and 1%) on the simulated price, pricePrediction  of the underlying security.


References

[2] Black-Scholes Model: What It Is, How It Works
[3] Scala for the Impatient - C. Horstman - Addison-Wesley 2012  

Friday, January 31, 2014

Performance Evaluation of Scala Maps

Target audience: Beginner

Overview

The Scala programming language API provides developers with a set of short immutable maps of predefined length, Map1, .. Map4,  Set. I thought it would be worthwhile to evaluate the performance of those short,  dedicated collections compare to their generic counterpart.
A second test consists of comparing HashMap with TreeMap for sorting a hash table by key.

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


Benchmark short maps

The benchmark consists of a very simple methods calls on immutable Map and Map4 instances as well as a mutable HashMap as defined in the code snippet below.

val map4 = new Map4("0"-> 0, "1" ->1, "2 ->,2, "3"-> 3)

val map0 = Map[String, Int]("0" -> 0, "1" -> 1, "2" -> 2, "3"-> 3)

val hMap = HashMap[String, Int]("0" -> 0, "1" -> 1, "2" -> 2, "3" -> 3)


We evaluate the time needed to execute 10 million iterations of a map, get and sum methods as follows:

aMap map { kv => kv._2 + 2 }
aMap get("2")
aMap values sum

The results of the performance test is shown in the table below. As expected the "short" collection is faster that the immutable Map and the mutable HashMap. However the performance improvement is not very significant.

Methods immutable.Map.Map4 immutable.Map mutable.HashMap
get 0.835 s 0.879 s 1.091 s
map 7.462 s 7.566 s 1.106 s
foldLeft 3.444 s 3.621 s 4.782 s


Sorting tables: TreeMap vs. sortBy

The second test consists of sorting a very large list or array of tuple (String, Float) by using maps. There are few options to sort a table among them:
  • Creating, populating and sorting a scala.collection.mutable.HashMap
  • Creating and populating a scala.collection.immutable.TreeMap
The test is set-up by creating pseudo-random key by concatenating the name of a city with an unique id. The values in the map are completely random.

 // Hashmap to be sorted 
val hMap = Range(0, sz)./:(new mutable.HashMap[String, Float])(
  (h, n) =>
      h += (s"${cities(Random.nextInt(cities.size))}_$n", Random.nextFloat )
)
val sortedByKey = hMap.toSeq.sortBy(_._1)

   // TreeMap
val treeMap = Range(0, sz)./:(new immutable.TreeMap[String, Float])(
    (h, n) => 
       h + ((s"${cities(Random.nextInt(cities.size))}_$n",  Random.nextFloat))
)

val sorted = sortedByKey.toSeq


The test is run with map which size varies between 10,000,00 and 90,000,00 entries


Sorting a tuple (String, Float) using TreeMap is roughly 40% faster than populating a HashMap and sort by key.

References

Monday, January 6, 2014

Scala Zip vs. Zipped

Target audience: Beginner
Estimated reading time: 3'

Overview

It is not unusual that Scala developers struggle in re-conciliating elegant functional programming style with efficient and fast execution. High order collection methods are conducive to very expressive constructs at the cost of poor performance. the zip method is no exception.
   def GenIterable.zip[B](that: GenIterable[B]): CC[(A, B)]
Fortunately, the authors of the Scala library have been diligent enough to provide us with an alternative in the case of the array of pairs (type Tuple2).

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

scala.Tuple2.zipped

Contrary to the GenIterable.zip method, Tuple2.zipped is a method unique and customized to the class Tuple2.

Let's evaluate the performance of the zipped relative to the ubiquitous zip. To this purpose, let's create a benchmark class to access elements on a zipped array. The first step is to create a function to access the first and last element for zipped arrays.
The first method zip exercises the GenIterable.zip method (line 2). The second method tzipped wraps the Tuple2.zip method (line 7).

1
2
3
4
5
6
7
8
9
def tzip(x: Array[Double], y: Array[Double]): (Double, Double) = {
    val z = x.zip(y)
   (z.head._1, z.last._2)
}
  
def tzipped(x: Array[Double], y: Array[Double]): (Double, Double) = {
    val z = (x, y).zipped
    (z.head._1, z.last._2)
}

Next we need to create a method that executes the two wrappers _zip and _zipped and measure the duration of their execution. We arbitrary select to sum the first element and the product of the last element of the zipped array.
The function zipTest has three arguments
  • A dataset of type Array[Double] x (line 2)
  • A second dataset y (line 3)
  • A function argument f for the wrapper methods tzip and tzipped (line 4)
The computation consists of computing the sum of elements for the first element of the tuple (line 12) and the product of the second element in the tuple (line 13).


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
def zipTest(
  x: Array[Double], 
  y: Array[Double], 
  f: (Array[Double], Array[Double]) => (Double, Double)
): Unit = {
  var startTime = System.currentTimeMillis
  
  var sum = 0.0
  var prod = 1.0
  Range(0, 50).foreach( _ => {
      val res = f(x, y)
      sum += res._1
      prod *= res._2
  })

  println(s"sum=$sum prod=$prod in ${(System.currentTimeMillis - startTime)}")


The last step is to invoke the benchmark method zipTest for different length of arrays. The element of the arrays are random floating point values..

1
2 
3
4
5
6
def zipTest(len: Int): Unit = {
    val x =  Array.tabulate(len)(_ => Random.nextDouble)
    val y = x.clone
    zipTest(x, y, _zip)
    zipTest(x, y, _zipped)
}


Performance Results

The test is performed on a single 8-core i7 32-Gbyte host.

val step = 1000000
val initial = step
Range(0, 6).foreach(n => zipTest(initial + n*step) )


















The results shows that the performance of Array.zip decreases exponentially compared to Tuple2.zipped which has a linear degradation. For 50 iterations on 1 million element arrays, Tuple2.zipped is 17% faster than Array.zip but 280% faster for 8 million elements.

Wednesday, December 18, 2013

Reinforcement Learning in Scala: Model

Target audience: Advanced
Estimated reading time: 5'

This post is the second installment of our introduction to reinforcement learning using the temporal difference. This section is dedicated to the ubiquitous Q-learning algorithm.

Table of contents
Follow me on LinkedIn

Overview

The previous post, Reinforcement learning in Scala: Policies introduced the concept of reinforcement learning, temporal difference and more specifically the Q-learning algorithm:
Reinforcement Learning I: States & Policy The following components were implemented:
  • States QLState and their associated actions QLAction
  • States space QLSpace
  • Policy QLPolicy as a container of tuple {reward, Q-value, probabilities} of the Q-learning model
The last two components to complete the implementation of Q-learning are the model and the training algorithm.

Modeling and training

The first step is to define a model for the reinforcement learning. A model is created during training and is composed of
  • Best policy to transition from any initial state to a goal state
  • Coverage ratio as defined as the percentage of training cycles that reach the (or one of the) goal.
class QLModel[T](val bestPolicy: QLPolicy[T], val coverage: Double)


The QLearning class takes 3 arguments
  • A set of configuration parameters config
  • The search/states space qlSpace
  • The initial policy associated with the states (reward and probabilities) qlPolicy

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
class QLearning[T](
   config: QLConfig, 
   qlSpace: QLSpace[T], 
   qlPolicy: QLPolicy[T]) 

    //model in Q-learning algorithm
  val model: Option[QLModel[T]] = train.toOption
    
    // Generate a model through multi-epoch training
  def train: Try[Option[QLModel[T]]] {}
  private def train(r: Random): Boolean {}

   // Predict a state as a destination of this current 
   // state, given a model
  def predict : PartialFunction[QLState[T], QLState[T]] {}

  // Select next state and action index
  def nextState(st: (QLState[T], Int)): (QLState[T], Int) {} 
}

The model of type Option[QLModel] (line 7) is created by the method train (line 10). Its value is None if training failed.

The training method train consists of executing config.numEpisodes cycle or episode of a sequence of state transition (line 5). The random generator r is used in the initialization of the search space.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def train: Option[QLModel[T]] = {
  val r = new Random(System.currentTimeMillis)

  Try {
    val completions = Range(0, config.numEpisodes).filter(train(r) )

    val coverage = completions.toSize.toDouble/config.numEpisodes
    if(coverage > config.minCoverage) 
       new QLModel[T](qlPolicy, coverage)
    else 
       QLModel.empty[T]
  }.toOption
}

The training process exits with the model if the minimum minCoverage (number of episodes for which the goal state is reached) is met (line 8).

The method train(r: scala.util.Random) uses a tail recursion to transition from the initial random state to one of the goal state. The tail recursion is implemented by the search method (line 4). The method implements the recursive temporal difference formula introduced in the previous post (Reinforcement Learning I: States & Policy) (lines 14-18). 

The state for which the action generates the highest reward R given a policy qlPolicy (line 10) is computed for each new state transition. The Q-value of the current policy is then updated qlPolicy.setQ before repeating the process for the next state, through recursion (line 21).

 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
30
31
32
33
def train(r: Random): Boolean = {
   
  @scala.annotation.tailrec
  def search(st: (QLState[T], Int)): (QLState[T], Int) = {
    val states = qlSpace.nextStates(st._1)
    if( states.isEmpty || st._2 >= config.episodeLength ) 
        (st._1, -1)
    
    else {
      val state = states.maxBy(s => qlPolicy.R(st._1.id, s.id))
      if( qlSpace.isGoal(state) )
          (state, st._2)
      else {
        val r = qlPolicy.R(st._1.id, state.id)   
        val q = qlPolicy.Q(st._1.id, state.id)
        // Q-Learning formula
        val deltaQ = r + config.gamma*qlSpace.maxQ(state, qlPolicy) -q
        val nq = q + config.alpha*deltaQ
        
        qlPolicy.setQ(st._1.id, state.id,  nq)
        search((state, st._2+1))
       }
     }
  } 
   
  r.setSeed(System.currentTimeMillis*Random.nextInt)

  val finalState = search((qlSpace.init(r), 0))
  if( finalState._2 != -1) 
     qlSpace.isGoal(finalState._1) 
  else 
     false
}

Note: There is no guaranty that one of the goal state is reached from any initial state chosen randomly. It is expected that some of the training epoch fails. This is the reason why monitoring coverage is critical. Obviously, you may choose a deterministic approach to the initialization of each training epoch by picking up any state beside the goal state(s), as a starting state.

Prediction

Once trained, the model is used to predict the next state with the highest value (or probability) given an existing state. The prediction is implemented as a partial function.

1
2
3
4
def predict : PartialFunction[QLState[T], QLState[T]] = {
  case state: QLState[T] if(model != None) => 
    if( state.isGoal) state else nextState(state, 0)._1
}

The method nextState does the heavy lifting. It retrieves the list of states associated with the current state st through its actions set (line 2). The next most rewarding state qState is computed using the reward matrix R of the best policy of the QLearning model (lines 6 - 8).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def nextState(st: (QLState[T], Int)): (QLState[T], Int) =  {
  val states = qlSpace.nextStates(st._1)
  if( states.isEmpty || st._2 >= config.episodeLength) 
    st
  else {
    val qState = states.maxBy(
     s => model.map(_.bestPolicy.R(st._1.id, s.id))
           .getOrElse(-1.0)
    )

    nextState( (qState, st._2+1))
  }
}


References