Parallel Monto-Carlo options pricing¶
This notebook shows how to use ipyparallel
to do Monte-Carlo options pricing in parallel. We will compute the price of a large number of options for different strike prices and volatilities.
Problem setup¶
%matplotlib inline
import matplotlib.pyplot as plt
import sys
import time
import numpy as np
Here are the basic parameters for our computation.
price = 100.0 # Initial price
rate = 0.05 # Interest rate
days = 260 # Days to expiration
paths = 10000 # Number of MC paths
n_strikes = 6 # Number of strike values
min_strike = 90.0 # Min strike price
max_strike = 110.0 # Max strike price
n_sigmas = 5 # Number of volatility values
min_sigma = 0.1 # Min volatility
max_sigma = 0.4 # Max volatility
strike_vals = np.linspace(min_strike, max_strike, n_strikes)
sigma_vals = np.linspace(min_sigma, max_sigma, n_sigmas)
from __future__ import print_function # legacy Python support
print("Strike prices: ", strike_vals)
print("Volatilities: ", sigma_vals)
Strike prices: [ 90. 94. 98. 102. 106. 110.]
Volatilities: [ 0.1 0.175 0.25 0.325 0.4 ]
Monte-Carlo option pricing function¶
The following function computes the price of a single option. It returns the call and put prices for both European and Asian style options.
def price_option(S=100.0, K=100.0, sigma=0.25, r=0.05, days=260, paths=10000):
"""
Price European and Asian options using a Monte Carlo method.
Parameters
----------
S : float
The initial price of the stock.
K : float
The strike price of the option.
sigma : float
The volatility of the stock.
r : float
The risk free interest rate.
days : int
The number of days until the option expires.
paths : int
The number of Monte Carlo paths used to price the option.
Returns
-------
A tuple of (E. call, E. put, A. call, A. put) option prices.
"""
import numpy as np
from math import exp,sqrt
h = 1.0/days
const1 = exp((r-0.5*sigma**2)*h)
const2 = sigma*sqrt(h)
stock_price = S*np.ones(paths, dtype='float64')
stock_price_sum = np.zeros(paths, dtype='float64')
for j in range(days):
growth_factor = const1*np.exp(const2*np.random.standard_normal(paths))
stock_price = stock_price*growth_factor
stock_price_sum = stock_price_sum + stock_price
stock_price_avg = stock_price_sum/days
zeros = np.zeros(paths, dtype='float64')
r_factor = exp(-r*h*days)
euro_put = r_factor*np.mean(np.maximum(zeros, K-stock_price))
asian_put = r_factor*np.mean(np.maximum(zeros, K-stock_price_avg))
euro_call = r_factor*np.mean(np.maximum(zeros, stock_price-K))
asian_call = r_factor*np.mean(np.maximum(zeros, stock_price_avg-K))
return (euro_call, euro_put, asian_call, asian_put)
We can time a single call of this function using the %timeit
magic:
%time result = price_option(S=100.0, K=100.0, sigma=0.25, r=0.05, days=260, paths=10000)
result
CPU times: user 174 ms, sys: 5.96 ms, total: 180 ms
Wall time: 211 ms
(12.166236181100073,
7.6440122060909745,
6.8409606562778666,
4.4639448434953621)
Parallel computation across strike prices and volatilities¶
The Client is used to setup the calculation and works with all engines.
import ipyparallel as ipp
rc = ipp.Client()
A LoadBalancedView
is an interface to the engines that provides dynamic load
balancing at the expense of not knowing which engine will execute the code.
view = rc.load_balanced_view()
Submit tasks for each (strike, sigma) pair. Again, we use the %%timeit
magic to time the entire computation.
async_results = []
%%time
for strike in strike_vals:
for sigma in sigma_vals:
# This line submits the tasks for parallel computation.
ar = view.apply_async(price_option, price, strike, sigma, rate, days, paths)
async_results.append(ar)
rc.wait(async_results) # Wait until all tasks are done.
CPU times: user 127 ms, sys: 15.7 ms, total: 143 ms
Wall time: 1.97 s
len(async_results)
30
Process and visualize results¶
Retrieve the results using the get
method:
results = [ar.get() for ar in async_results]
Assemble the result into a structured NumPy array.
prices = np.empty(n_strikes*n_sigmas,
dtype=[('ecall',float),('eput',float),('acall',float),('aput',float)]
)
for i, price in enumerate(results):
prices[i] = tuple(price)
prices.shape = (n_strikes, n_sigmas)
Plot the value of the European call in (volatility, strike) space.
plt.figure()
plt.contourf(sigma_vals, strike_vals, prices['ecall'])
plt.axis('tight')
plt.colorbar()
plt.title('European Call')
plt.xlabel("Volatility")
plt.ylabel("Strike Price")
<matplotlib.text.Text at 0x106036470>

Plot the value of the Asian call in (volatility, strike) space.
plt.figure()
plt.contourf(sigma_vals, strike_vals, prices['acall'])
plt.axis('tight')
plt.colorbar()
plt.title("Asian Call")
plt.xlabel("Volatility")
plt.ylabel("Strike Price")
<matplotlib.text.Text at 0x1062d27f0>

Plot the value of the European put in (volatility, strike) space.
plt.figure()
plt.contourf(sigma_vals, strike_vals, prices['eput'])
plt.axis('tight')
plt.colorbar()
plt.title("European Put")
plt.xlabel("Volatility")
plt.ylabel("Strike Price")
<matplotlib.text.Text at 0x106c0c748>

Plot the value of the Asian put in (volatility, strike) space.
plt.figure()
plt.contourf(sigma_vals, strike_vals, prices['aput'])
plt.axis('tight')
plt.colorbar()
plt.title("Asian Put")
plt.xlabel("Volatility")
plt.ylabel("Strike Price")
<matplotlib.text.Text at 0x1078690b8>
