# Parallel Monte Carlo with Dask

Today, we will look at how to implement a Monte Carlo simulation
to value an Asian option in a Jupyter notebook; and how to distribute that Monte Carlo
simluation across a cluster using Dask. Don’t worry if you are not a
Quantitative Analyst, this blog is more about Dask, `from dask import delayed`

and exciting
cluster stuff than about financial mathematics.

## Asian options

An Asian option is an exotic type of financial option contract. The payoff at expiry
depends on the average underlying price over some pre-set period of time - different to
a plain vanilla option where the payoff depends on the underlying price at expiry (fun
fact: they are called *Asian* options because legend has it that
Bankers Trust traded the first option linked to the average price of crude oil in in Tokyo
back in 1987).

The payoff of a fixed-strike Asian call option at expiry is defined as: `C(T) = max(0, A(T) - K)`

.
There is no closed form analytical formula to determine the fair value any time before
expiry for this kind of option. But we can determine the fair value by simulating a
large number of paths for the underlying asset price and calculate the payout at expiry
for each path. The fair value is then the discounted average over all these Monte Carlo
simulations.

## Random walk

We assume that the current value of the underlying asset price is composed of the past value and
a determinstic upward trend plus an error term defined as a white noise (a normal variable with zero mean and variance one). This is called a random walk with drift. Given a start price `s0`

, a constant drift `mu`

and volatility `sigma`

, we can simulate the price path over `days`

number of days like this:

```
import numpy as np
def random_walk(s0, mu, sigma, days):
dt = 1/365.
prices = np.zeros(days)
shocks = np.zeros(days)
prices[0] = s0
for i in range(1, days):
e = np.random.normal(loc=mu * dt, scale=sigma * np.sqrt(dt))
prices[i] = prices[i-1] * (1 + e)
return prices
```

For any random path, we can now calculate the average price and, given the strike price `K`

, the payoff at expiry:

```
days = 365 * 4 # days to expiry
s0 = 100 # current underlying price
mu = 0.02 # drift
sigma = 0.2 # volatility
K = 100 # strike price
A = np.average(random_walk(s0, mu, sigma, days)
C = max(0, A - K)
```

## Monte Carlo Simulation

We need to generate a large number of random price paths for the underlying. For each price path we calculate the associated payoff. These payoffs are averaged and discounted to today. This result is the value of the option.
The actual number of simulations `n`

required depends on the confidence level you are comfortable with for
the estimate. That, in turn, depends on your option’s parameters. I won’t go into details here, so let’s just assume that `n = 10000`

. We need to generate `n`

random paths and calculate the average value of the option at expiry:

```
n = 10000
%time np.average([max(0, np.average(random_walk(s0, mu, sigma, days)) - K) for i in range(0, n)])
Wall time: 1min 10s
11.673350929139112
```

## Dask task scheduler: dask.distributed

If you have more than one CPU at your disposal, you can bring down the calculation time by
distributing the random walk generation across multiple CPUs. This is where Dask comes in.
Dask is a flexible library for parallel computing in Python, optimised for
interactive computational workloads. Instead of running `n`

random walks in a single
Jupyter notebook process, we make Dask distribute these calculations across
a number of processes.

Distributing these calculations is the job of the `task scheduler`

. Think of the task
scheduler as a server that manages where and how calculations as executed. There are different types of schedulers: single machine scheduler (limited to local machine) and distributed scheduler (distributed across a cluster). Here, I will
use a simple single machine scheduler.

In a new shell/command window, start a local `dask.distributed`

process scheduler:

```
~$ python
...
>>> from dask.distributed import Client
>>> client = Client(scheduler_port=8786)
```

You can interrogate the `client`

object to confirm it runs on port 8786 on localhost.
By default, it spins up as many processes as you have CPUs on your computer, in my case I have 4 cores:

```
>>> client
<Client: scheduler='tcp://127.0.0.1:8786' processes=4 cores=4>
```

## Wrapping functions with dask.delayed

We need to instruct our `random_walk`

and `np.average`

function calls
to execute lazily. Instead of executing the functions immediately,
we want to defer execution via the Dask task scheduler. We can achieve
this by wrapping functions in `dask.delayed`

. Alternatively, you can decorate
your functions with `@dask.delayed`

(not covered in this blog post, though
have a look at http://dask.pydata.org/en/latest/delayed.html#decorator).
This delays the execution of the function and generates a task graph instead.

```
from dask import delayed
s0 = 100
K = 100
mu = 0.02
sigma = 0.2
days = 365*4
n = 10000
result = delayed(np.average)([
delayed(max)(
0,
delayed(np.average)(random_walk(s0, mu, sigma, days)) - K
) for i in range(0, n)
])
```

`result`

is now a `Delayed`

object and contains the task graph. The task graph is a
directed acyclic graph (DAG) and models the dependencies between the `np.average`

,
`max`

and `random_walk`

function calls - without executing them. This allows Dask to optimise
its task execution strategy.

## Execute the task graph

We have everything in place now to execute the actual calculation of `result`

. Remember
that `result`

up to this point is just a Delayed object. Execute the computation:

```
%time result.compute()
Wall time: 39.9 s
11.617145839123664
```

In essence, the calculation time has roughly halved on 4 cores compared to 1 core. The reason it has not reduced to 25% is that there is some overhead involved in managing and distributing the task themselves. Have a look at the Dask docs about dask.delayed and task scheduling.