An Introduction to Monte Carlo: A Useful Tool to Estimate Integrals and Option Prices

Often in numerical analysis, it can be very complex to find exact solutions to problems. Fortunately, exact answers are often not needed and estimations can be used to make a problem much simpler to get a close answer. One such technique is through Monte Carlo, which relies on repeatedly sampling random values from a probability distribution to perform an estimation. With enough samples, the estimated values converge to a mean based on the law-of-large-numbers.

Law Of Large Numbers Definition

A good use case of Monte Carol is estimating integrals. Recall that an expected value of a random variable E(x) which has a probability density function f(x) can be defined as:

The integral bounds depend on the range of the probability function

Thus, if we are given an integral which we want to estimate, we can easily convert it to an expectation instead of an integral. For example, if we want to evaluate to integrate sin(x) from 0 to π, we can rewrite this calculation as:

The 1/π is the probability density function and is a uniform distribution

We can set up Monte Carlo Estimation through repeatedly sampling an uniform distribution. If we input samples uniformally from 0 to π into sin(x) and take the mean of the values, that will give us an estimate for the integral. Monte Carlo sounds complicated but really, it comes down to sampling from a distribution, plugging into a equation repeatedly and taking the average. Here is a code snippet that shows the calculation here:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import uniform
msin_estimates = [None]*99
msin_std = [None] *99
for i in range(1,100):
    #sample uniformly numbers from 0 to pi
    unif_array = uniform.rvs(size = i*500)*np.pi
    # Plug into equation above
    sin_val = np.sin(unif_array)*np.pi
    # Calculate the mean which we use as the estimate
    msin_estimates[i-1] = np.mean(sin_val)
    msin_std[i-1] = np.std(sin_val)/np.sqrt(i*500)

Now, let’s graph it. We can plot these values against the actual value of the integral, 2. Also, let’s plot the standard deviation as well to show that as we increase the sample size, the standard deviation converges to 0 as the mean converges to the actual value.

## Graph Generation
import matplotlib.pyplot as plt
plt.plot([2]*99)
plt.plot(msin_estimates, '.')
plt.plot(2+np.array(msin_std)*3)
plt.plot(2+np.array(msin_std)*-3)
plt.xlabel('Sample size, per 500')
plt.ylabel('Value')

Another application of Monte Carlo is option pricing in finance where random walks of a share price are simulated and used to determine how much an option is worth using the Black-Scholes equation for prices for call and put options, whose formula can be found here:

Black-Scholes Formula (d1, d2, Call Price, Put Price, Greeks)

Let’s estimate a price for a European put option with the following information

• Share volatility is 30%
• The risk-free rate is 10%
• Average share return is 15%
• The current share price is $100
• The strike on the European put option is $110
• The term of the contract is 6 months

First, we use two helper functions: terminal_value_of_share simulates a random walk and we are simulating Z, which is from the normal distribution and gives the z-score, meaning how many standard deviations away an estimate is from its mean. From that we can plug into the equation that gives us the price of discounted put payoff, which factors the risk_free_rate as well.

def terminal_value_of_share(risk_free_rate, S_0, sigma,Z, T):
    return S_0*np.exp((risk_free_rate-sigma**2/2)*T+sigma*np.sqrt(T)*Z)
def discounted_put_payoff(S_T, K, risk_free_rate, T):
    return np.exp(-risk_free_rate*T)*np.maximum(K-S_T,0)

For this estimation, we a sampling a normal distribution and estimated what the terminal value of share would and from there getting the price of the put contract. I have also calculated the analytical price directly from the Black-Scholes equation to see how they compare.

from scipy.stats import norm
import math
risk_free =.1
S_0 = 100
sigma = .3

strike = 100
T = .5
current_time = 0
put_estimates = [None] *50
put_std = [None]*50
for i in range(1,51):
    norm_array = norm.rvs(size=1000*i)
    term_val= terminal_value_of_share( risk_free, S_0 ,sigma, norm_array, T)
    put_val = discounted_put_payoff(term_val, strike, risk_free, T)
    put_estimates[i-1] = np.mean(put_val)
    put_std[i-1] = np.std(put_val)/np.sqrt(i*1000)
d_1 = (math.log(S_0/strike) + (risk_free +sigma**2/2)* (T))/(sigma*math.sqrt(T))
d_2 = d_1 -sigma*math.sqrt(T)
analytic_putprice = -S_0*norm.cdf(-d_1) + strike*math.exp(-risk_free*(T-current_time))*norm.cdf(-d_2)

Plotting the put prices estimates against the theoretical value:

plt.plot(put_estimates, '.')
plt.plot([analytic_putprice]*50)
plt.plot(analytic_putprice + 3 *np.array(put_std), 'r')
plt.plot(analytic_putprice -3 *np.array(put_std), 'r')
plt.xlabel("Sample Size")
plt.ylabel("Value")
plt.show()

As you can see, the estimates again converge to the theoretical value. The examples I have shown are straightforward to calculate but Monte Carlo becomes indispensable when dealing with dealing with problems where it is difficult to calculate the exact solution.