Permalink
Browse files

update parallel demos for newparallel

  • Loading branch information...
minrk committed Feb 16, 2011
1 parent 5bcc966 commit afabc88a82ebde75995e044bb0263d8ea18cac95
@@ -0,0 +1,144 @@
+#!/usr/bin/env python
+"""Run a Monte-Carlo options pricer in parallel."""
+
+#-----------------------------------------------------------------------------
+# Imports
+#-----------------------------------------------------------------------------
+
+import sys
+import time
+from IPython.zmq.parallel import client
+import numpy as np
+from mcpricer import price_options
+from matplotlib import pyplot as plt
+
+#-----------------------------------------------------------------------------
+# Setup parameters for the run
+#-----------------------------------------------------------------------------
+
+def ask_question(text, the_type, default):
+ s = '%s [%r]: ' % (text, the_type(default))
+ result = raw_input(s)
+ if result:
+ return the_type(result)
+ else:
+ return the_type(default)
+
+cluster_profile = ask_question("Cluster profile", str, "default")
+price = ask_question("Initial price", float, 100.0)
+rate = ask_question("Interest rate", float, 0.05)
+days = ask_question("Days to expiration", int, 260)
+paths = ask_question("Number of MC paths", int, 10000)
+n_strikes = ask_question("Number of strike values", int, 5)
+min_strike = ask_question("Min strike price", float, 90.0)
+max_strike = ask_question("Max strike price", float, 110.0)
+n_sigmas = ask_question("Number of volatility values", int, 5)
+min_sigma = ask_question("Min volatility", float, 0.1)
+max_sigma = ask_question("Max volatility", float, 0.4)
+
+strike_vals = np.linspace(min_strike, max_strike, n_strikes)
+sigma_vals = np.linspace(min_sigma, max_sigma, n_sigmas)
+
+#-----------------------------------------------------------------------------
+# Setup for parallel calculation
+#-----------------------------------------------------------------------------
+
+# The Client is used to setup the calculation and works with all
+# engines.
+c = client.Client(profile=cluster_profile)
+
+# 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 = c[None]
+
+# Initialize the common code on the engines. This Python module has the
+# price_options function that prices the options.
+
+#-----------------------------------------------------------------------------
+# Perform parallel calculation
+#-----------------------------------------------------------------------------
+
+print "Running parallel calculation over strike prices and volatilities..."
+print "Strike prices: ", strike_vals
+print "Volatilities: ", sigma_vals
+sys.stdout.flush()
+
+# Submit tasks to the TaskClient for each (strike, sigma) pair as a MapTask.
+t1 = time.time()
+async_results = []
+for strike in strike_vals:
+ for sigma in sigma_vals:
+ ar = view.apply_async(price_options, price, strike, sigma, rate, days, paths)
+ async_results.append(ar)
+
+print "Submitted tasks: ", len(async_results)
+sys.stdout.flush()
+
+# Block until all tasks are completed.
+c.barrier(async_results)
+t2 = time.time()
+t = t2-t1
+
+print "Parallel calculation completed, time = %s s" % t
+print "Collecting results..."
+
+# Get the results using TaskClient.get_task_result.
+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)
+strike_mesh, sigma_mesh = np.meshgrid(strike_vals, sigma_vals)
+
+print "Results are available: strike_mesh, sigma_mesh, prices"
+print "To plot results type 'plot_options(sigma_mesh, strike_mesh, prices)'"
+
+#-----------------------------------------------------------------------------
+# Utilities
+#-----------------------------------------------------------------------------
+
+def plot_options(sigma_mesh, strike_mesh, prices):
+ """
+ Make a contour plot of the option price in (sigma, strike) space.
+ """
+ plt.figure(1)
+
+ plt.subplot(221)
+ plt.contourf(sigma_mesh, strike_mesh, prices['ecall'])
+ plt.axis('tight')
+ plt.colorbar()
+ plt.title('European Call')
+ plt.ylabel("Strike Price")
+
+ plt.subplot(222)
+ plt.contourf(sigma_mesh, strike_mesh, prices['acall'])
+ plt.axis('tight')
+ plt.colorbar()
+ plt.title("Asian Call")
+
+ plt.subplot(223)
+ plt.contourf(sigma_mesh, strike_mesh, prices['eput'])
+ plt.axis('tight')
+ plt.colorbar()
+ plt.title("European Put")
+ plt.xlabel("Volatility")
+ plt.ylabel("Strike Price")
+
+ plt.subplot(224)
+ plt.contourf(sigma_mesh, strike_mesh, prices['aput'])
+ plt.axis('tight')
+ plt.colorbar()
+ plt.title("Asian Put")
+ plt.xlabel("Volatility")
+
+
+
+
+
+
@@ -0,0 +1,45 @@
+
+def price_options(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)
+
@@ -0,0 +1,63 @@
+"""Calculate statistics on the digits of pi in parallel.
+
+This program uses the functions in :file:`pidigits.py` to calculate
+the frequencies of 2 digit sequences in the digits of pi. The
+results are plotted using matplotlib.
+
+To run, text files from http://www.super-computing.org/
+must be installed in the working directory of the IPython engines.
+The actual filenames to be used can be set with the ``filestring``
+variable below.
+
+The dataset we have been using for this is the 200 million digit one here:
+ftp://pi.super-computing.org/.2/pi200m/
+
+and the files used will be downloaded if they are not in the working directory
+of the IPython engines.
+"""
+
+from IPython.zmq.parallel import client
+from matplotlib import pyplot as plt
+import numpy as np
+from pidigits import *
+from timeit import default_timer as clock
+
+# Files with digits of pi (10m digits each)
+filestring = 'pi200m.ascii.%(i)02dof20'
+files = [filestring % {'i':i} for i in range(1,16)]
+
+# Connect to the IPython cluster
+c = client.Client()
+c.run('pidigits.py')
+
+# the number of engines
+n = len(c.ids)
+id0 = list(c.ids)[0]
+# fetch the pi-files
+print "downloading %i files of pi"%n
+c.map(fetch_pi_file, files[:n])
+print "done"
+
+# Run 10m digits on 1 engine
+t1 = clock()
+freqs10m = c[id0].apply_sync_bound(compute_two_digit_freqs, files[0])
+t2 = clock()
+digits_per_second1 = 10.0e6/(t2-t1)
+print "Digits per second (1 core, 10m digits): ", digits_per_second1
+
+
+# Run n*10m digits on all engines
+t1 = clock()
+c.block=True
+freqs_all = c.map(compute_two_digit_freqs, files[:n])
+freqs150m = reduce_freqs(freqs_all)
+t2 = clock()
+digits_per_second8 = n*10.0e6/(t2-t1)
+print "Digits per second (%i engines, %i0m digits): "%(n,n), digits_per_second8
+
+print "Speedup: ", digits_per_second8/digits_per_second1
+
+plot_two_digit_freqs(freqs150m)
+plt.title("2 digit sequences in %i0m digits of pi"%n)
+plt.show()
+
Oops, something went wrong.

0 comments on commit afabc88

Please sign in to comment.