Monte Carlo Techniques: Estimating the Value of Pi with Random Numbers

I want to write a series of articles on using Monte Carlo methods in computational finance.  Monte Carlo algorithms are numerical techniques that use repeated sampling of random numbers to perform the desired calculation.  Such techniques are becoming increasingly frequent even down to the level of the individual trader with some brokerages using Monte Carlo simulations to estimate the probability of reaching a given percentage of maximum profit on short-premium options trades.

While that is not as complicated as one might think, I wanted to start out with a very simple article first.  In this article, we will write a simple Python script to estimate the value of $\pi$. In a follow-up post, we will introduce the concept of a random walk and use it to model a stock's price action. In the last article of this series, we will show how to predict the probability of reaching a given percentage of profit on an option trade.

So how can we use random numbers to estimate a constant like $\pi$?  We start with a unit circle circumscribed by a square.  Since the circle has a radius of one, the side of the square must have a length of two.

Figure 1:  A unit circle centered at the origin circumscribed by a square.

For the sake of coding convinces, we'll only look at the upper right quadrant as shown below in Fig 2.  This is because we'll use the built in rand function in numpy, which selects random values between zero and one. If we restrict ourselves in this way, we won't have to rescale this range to between -1 and 1.  Rescaling is no big deal, but if we can simply avoid doing so, we will take that opportunity.

Figure 2:  We will only consider the upper left quadrant.

The square shaded above in light gray now has sides $L = 1$ and an area of $L^2 = 1$.  The quadrant in the shaded region has an area of one-quarter the original, in other words, $\pi R^2 / 4$.  If we throw a large number of darts that hit randomly in that upper right quadrant,  we would expect the ratio of the number that hit within the circle to the total number of darts thrown to approach the ratio of the circle's area to that of the square.  In other words,
$$\frac{N_{inside}}{N_{total}} \approx \frac{A_{circle}}{A_{square}} = \frac{ \frac{1}{4} \pi R^2}{1}.$$

Fig. 3 shows an animation of this process.  For the sake of clarity, we only include 25 points in this animation whereas the actual calculation we will use significantly more.

Figure 3:  25 XY-pairs drawn at random.  For large numbers of pairs, the ratio of those within the circle to the total number of points should converge to the ration of the circle's area to that of the square.

Finding if a given point is within the circle is straightforward.  Since the circle has a radius of 1, the point is inside if
$$\sqrt{ x^2 + y^2} \leq 1^2$$.
The code below implements this calculation using the Python's Numpy library.  It should print "pi =  3.140936" to the screen.

# Python 2.7 script to estimate the value of pi using Monte Carlo techniques

import numpy as np

#  Use a constant random seed for debugging purposes
np.random.seed(0)

#  Number of random points to choose
N = 1000000

#  Create a 2-dimensional numpy array.  Each row will correspond to one
#  point with the entry in the first column being the X value and the
#  second column providing the Y values.
p = np.random.rand(N, 2)

#  Use the built-in norm function to calculate the distance from the origin
n = np.linalg.norm(p, axis = 1, keepdims = True)

#  Sum up the number of entries in the circle.  That is count the number of
#  points who have a norm less than or equal to one
number_in_circle =  np.sum(n <= 1)

#  Print out the results
print 'pi = ', 4 * float(number_in_circle) / float(N)

While this example borders on the trivial, it serves as a good starting point for more sophisticated calculations.  As noted in the introduction, I've used Monte Carlo methods in financial modeling (which we will discuss in later posts), as well as in complex physical systems such as modeling the loss of active material from cycling batteries.

Some comments on vectorization

The code above uses a "vectorized" implementation.  If you are familiar with the concept, there is no need to read further.  If not, I'd like to take a moment to explain the concept because we will use it heavily in future articles.

Numpy is a package for numerical calculations in Python.  It provides a significant speed increase over writing code in straight Python because it is written in C rather than Python itself.  It also provides the ability to "vectorize" calculations meaning operation can be performed on whole vectors or matrices without having to explicitly loop through the elements.  This can further enhance speed considerably.

Figure 4:  Snapshot of a Jupyter notebook.  We load in the numpy and math libraries (the math functionality will be used later), create two numpy arrays, add them together, and finally print out the results.

As an example, in Fig 4, we create two arrays p1, and p2, and compute their elementwise sum.  There is no need to loop over each index and add them in a piecemeal manner.  Numpy also provides a number of functions that can perform vectorized operation on arrays.  Notice in the code to estimate $\pi$,  we use the functions np.linalg.norm and np.sum to calculate the distance of the point from the origin and count the number of points with the circle, respectively.

Figure 5:  Compute the L2 norm of a vector $p$, and compare the result to that of calculating the distance between $p$ and the point (0, 0) using the Pythagorean theorem.

The function np.linalg.norm calculates the length of a vector (L2 norm).  In the code snippet shown in Fig. 5 we show that calculating the L2 norm via this  functionality is the same as explicitly calculating the length between the origin and a point $p$.  Lastly in Fig 6. we show the equivalence of looping over an array and calculating the distance of each point in that array to simply using the built in norm function.

Figure 6:  Combining the notions from Figs 4 & 5, we compare the vectorized implementation of the norm of an array of points to that of explicitly looping over the array and calculating the L2 norm via the Pythagorean theorem.

The vectorization used here is straightforward and pretty intuitive, but it is important to have a handle on the concept.  In subsequent articles where we discuss ransom walks in the context of modeling stock behavior, we will lean heavily on this as well as SciPy's ability to deal with sparse matrices  rapidly.  This is a bit more abstract, and thus it is a good idea to make sure one understands the simple case before jumping into the deep end.

Articles in this series

Part II:  Monte Carlo Techniques: Modeling Stock Price Action
Part III:  Monte Carlo Techniques:  Calculating the Probability of Making 50% of Max Profit on Short Option Positions
Part IV:  Monte Carlo Techniques:  Estimating Correlations between Option Positions and  and Their Underlying Stock

Calculating Implied Volatility from an Option Price

Nowadays, with the ubiquity of computers and information, many traders wish to know how the sausage is made and work out for themselves some of the numbers appearing on the screen of their trading platforms.  This is often a simple matter of looking up the formulas involved and plugging in the numbers.

When calculating the numbers pertaining to options, however, we run into an issue with implied volatility.  The Black-Scholes model tells us what an option should be worth given its strike price, the risk-free interest rate, the remaining time until expiration, the stock's price, and the implied volatility.  For example, the price of a call, $C$, is given by,

$$C = \Phi(d_1) S - \Phi(d_2) K e^{-r t},$$

where $S$ is the price of the stock, $K$ is the strike price, $r$ is the annualized risk-free rate, and $t$ is the remaining time to expiration expressed in years.  $\Phi$ is the normal cumulative distribution function, and $d_1$ and $d_2$ are given by,
$$d_1 = \frac{1}{\sigma \sqrt{t}} \left[ \ln\left(\frac{S}{K}\right) + \left(r + \frac{\sigma.^2}{2}\right) t\right],$$
and
$$d_2 = d_1 - \sigma \sqrt{t},$$
respectively.  In the above expressions for $d_1$ and $d_2$,  $\sigma$ is the implied volatility.

If we wish to calculate $\sigma$, we run into an issue.  All the variables in the above equations are known, and we can get the call price directly from the option chain, but we are incapable of isolating $\sigma$ algebraically.  Instead, we will have to turn to numerical methods to calculate the implied volatility.

Recall from high school algebra that if $y = f(x)$, the value of $x$ for which $f(x) = 0$ is called the root of the function $f$.  We will make use of a root finding algorithm to find our volatility, $\sigma$.  Our function will be the theoretical call price from the Black-Scholes model minus the known option price.  We will insist that equals zero and find the value of $\sigma$ that makes it so.  Explicitly, we want,

$$C(\sigma) - C_0 = 0,$$

where $C_0$ is the call price from the option chain.

For the sake if this article, we will assume the underlying stock is trading for 100, the strike price is 105, there are 30 days until the contract expires, and the risk-free rate is 1%.  We will plot $C(\sigma) - C_0$ for values of implied volatility from zero to one below just to get a sense of how this function behaves.

Figure 1.  We wish to find the value of $\sigma$ where $C(\sigma) - C_0$ is zero.

We wish to find the value of $\sigma$ where the function $C(\sigma) - C_0 = 0$.  From Fig 1., we can see this happens around $\sigma = 0.38$.  We will use an iterative method for finding an approximation to the root developed by Isaac Newton and Joseph Raphson.  The method requires an initial guess to start the process.  We will denote our successive approximations of the root as $\sigma_i, i = 0,1,2,...$ and will start with a guess of $\sigma_0 = 0.5$.

The Newton-Raphson method calculates the slope of the tangent line evaluated at the first iteration point.  It then uses the point where that tangent line crosses zero as the starting point for the next iteration.  This is shown graphically in Fig. 2.

Figure 2.  Plot of the line tangent to our function at $\sigma = 0.5$.  The intercept of this line with the x-axis is used as the starting point for the next iteration.  We can see after only one iteration it is relatively close to the actual root

The process then repeats until our test value for $\sigma$ is as close to zero as we like.  Fig. 3 shows the second iteration.  Note that after only two iterations, we are very close to the root.

Figure 3.  Tangent line for the second iteration

To implement this technique, we need to be able to calculate the derivative of the function of interest.  One of the nice features of the Black-Scholes model is that it has a closed form that can be easily differentiated.  Traders refer to the derivative of the option price with respect to volatility as vega and denote it with the Greek letter $\nu$.  Vega is given by,
$$\nu = \frac{\partial C}{\partial\sigma} = S \phi(d_1)\sqrt{t},$$
where $\phi$ is the normal probability density function.  With the above equations, we have enough information to implement a program to calculate the implied volatility of an option.  We will use Python for this exercise because it is a popular, freely available programming language that has a fairly extensive math and statistics libraries.

We will make use of the scypi.stats library as well as specific functions in the math package.  After importing these, we need to write a couple of helper functions to implement the Black-Scholes model for call option prices.  The first function calculates the values of $d_1$ and $d_2$.  The second function actually calculates the theoretical price of the call.

from scipy.stats import norm
from math import sqrt, exp, log, pi

def d(sigma, S, K, r, t):
    d1 = 1 / (sigma * sqrt(t)) * ( log(S/K) + (r + sigma**2/2) * t)
    d2 = d1 - sigma * sqrt(t)
    return d1, d2

def call_price(sigma, S, K, r, t, d1, d2):
    C = norm.cdf(d1) * S - norm.cdf(d2) * K * exp(-r * t)
    return C

Next we enter in the known option information as well as our initial guess for $\sigma$.

#  S  = Stock price
#  K  = strike
#  C  = price of call as predicted by Black-Scholes model
#  r  = risk-free interest rate
#  t  = time to expiration expressed in years
#  C0 = price of call option from option chain

S = 100.0
K = 105.0
r = 0.01
t = 30.0/365
C0 = 2.30

#  We need a starting guess for the implied volatility.  We chose 0.5
#  arbitrarily.
vol = 0.5

Next we define some variables needed for bookkeeping.  These will include a variable to count the number of iterations along with a maximum number of iterations.  We will use these to make sure our program doesn't get stuck in an infinite loop.  We will also set a tolerance where we will terminate the iterations when our test root gives a value that is close enough to zero to satisfy our needs.

epsilon = 1.0          #  Define variable to check stopping conditions
abstol = 1e-4          #  Stop calculation when abs(epsilon) < this number

i = 0                  #  Variable to count number of iterations
max_iter = 1e3         #  Max number of iterations before aborting

Lastly, we will do the iteration and print out the results.

while epsilon > abstol:
    #  if-statement to avoid getting stuck in an infinite loop.
    if i > max_iter:
        print 'Program failed to find a root.  Exiting.'

        break

    i = i + 1
    orig = vol
    d1, d2 = d(vol, S, K, r, t)
    function_value = call_price(vol, S, K, r, t, d1, d2) - C0
    vega = S * norm.pdf(d1) * sqrt(t)
    vol = -function_value/vega + vol
    epsilon = abs(function_value)

print 'Implied volatility = ',  vol
print 'Code required', i, 'iterations.'

Running this produces the result,

Implied volatility =  0.368856324914
Code required 3 iterations.

We see the program converges quickly taking only three iterations to find a value within our tolerance.  The above technique can be used for puts by substituting in the Black-Scholes formula for put prices, of course.   So that's it.  Pretty simple.

Newton's method works well when the function and its derivative are well-behaved.  In our case, the function only has one root, so we don't have to worry about the possibility of having several mathematically allowable answers and setting up our code and initial guess to converge to the proper root or deciding which of several results is the right answer to our problem.

Here is the above code in its entirety.

#!/usr/bin/python

from scipy.stats import norm
from math import sqrt, exp, log, pi

def d(sigma, S, K, r, t):
    d1 = 1 / (sigma * sqrt(t)) * ( log(S/K) + (r + sigma**2/2) * t)
    d2 = d1 - sigma * sqrt(t)
    return d1, d2

def call_price(sigma, S, K, r, t, d1, d2):
    C = norm.cdf(d1) * S - norm.cdf(d2) * K * exp(-r * t)
    return C

#  S  = spot
#  K  = strike
#  C  = price of call as predicted by Black-Scholes model
#  r  = risk-free interest rate
#  t  = time to expiration expressed in years
#  C0 = price of call option from option chain

S = 100.0
K = 105.0
r = 0.01
t = 30.0/365
C0 = 2.30

#  We need a starting guess for the implied volatility.  We chose 0.5
#  arbitrarily.
vol = 0.5

epsilon = 1.0  #  Define variable to check stopping conditions
abstol = 1e-4  #  Stop calculation when abs(epsilon) < this number

i = 0   #  Variable to count number of iterations
max_iter = 1e3  #  Max number of iterations before aborting

while epsilon > abstol:
    #  if-statement to avoid getting stuck in an infinite loop.
    if i > max_iter:
        break

    i = i + 1
    orig = vol
    d1, d2 = d(vol, S, K, r, t)
    function_value = call_price(vol, S, K, r, t, d1, d2) - C0
    vega = S * norm.pdf(d1) * sqrt(t)
    vol = -function_value/vega + vol
    epsilon = abs(function_value)

print 'Implied volatility = ',  vol
print 'Code required', i, 'iterations.'

For further reading see:

Wikipedia article on the Black-Scholes model
Wikipedia article on Newton's method

Updated 11/12/2018:  Fixed text which said the strike price was 155 with a stock price of 150.  The code uses a strike of 105 and a price of 100.  Sorry for not catching this sooner.