European Vanilla Option Pricing - Monte Carlo Methods

We start with the assumption that underlying follow Geometric Brownian Motion (GBM):

$dS_{t} = \mu S_{t}dt + \sigma S_{t}dW_{t}$

We use Ito’s Lemma with $G = \ln (S)$ , then we have

$\frac{\partial G}{\partial S} = \frac{1}{S}$

$\frac{\partial ^{2}G}{\partial S^{2}} = -\frac{1}{S^{2}}$

$\frac{\partial G}{\partial t} = 0$

By Ito’s Lemma, we have

$dG = (\mu - \frac{\sigma ^{2}}{2})dt + \sigma dz$

Therefore, the change of $\ln(S)$ between time 0 and future time T, is normally distributed as following:

$\ln(S_{T}) - \ln(S_{0}) \sim \mathcal{N}[(\mu - \frac{\sigma ^{2}}{2})T, \sigma \sqrt{T}]$

Thus, the future underlying price can be written as,

$S_{T} = S_{0}e^{(\mu-\frac{\sigma ^{2}}{2})T + \sigma \sqrt{T}\epsilon}$

$\epsilon$ is the noise term from standard normal distribution.

Note, we will take $\mu = r$ , which is the risk free rate. This means investors are risk neutral and requires risk free return on underlying asset. This is to be consistent with the risk neutral probabilities used in simulation.
Correspondingly, we also use risk free rate in the discount factor in step 5.

So now we can simulate the future underlying price at expiry. With European Call or Put boundary condition to calculate the payoff.

$Call Payoff at Expiry = Max[(S_{T} - K), 0]$

$Put Payoff at Expiry = Max[(K - S_{T}), 0]$

We then need to discount the future payoff back to present by multiplying a discount factor,

$Discount Factor = e^{-rT}$

The above two steps are repeated many times and its expectation is calculated as the final simulation result.

python code：

def getMCPrice(self):
    """
    Determine the option price using a Monte Carlo approach.
    The log return of underlying follow Normal distribution.
    s_T = s_t * exp((r - 1/2 * sig^2) * (T-t) + sig * sqrt(T-t) * sig_Normal)
    """
    calc = np.zeros([self.iterations, 2])
    rand = np.random.normal(0, 1, [1, self.iterations])
    mult = self.spot * np.exp(self.tenor * (self.rate - 0.5 * self.sigma**2))
 
    if self.callPut == 'Call':
        calc[:,1] = mult * np.exp(np.sqrt((self.sigma**2)*self.tenor) * rand) - self.strike
    elif self.callPut == 'Put':
        calc[:,1] = self.strike - mult*np.exp(np.sqrt((self.sigma**2) * self.tenor) * rand)
 
    avgPayOff = np.sum(np.amax(calc, axis=1)) / float(self.iterations)
 
    return np.exp(-self.rate * self.tenor) * avgPayOff
 
def getBSPrice(self):
    """ Determine the option price using the exact Black-Scholes expression. """
    return blackScholesOptionPrice(self.callPut, self.spot, self.strike, self.tenor, self.rate, self.sigma)

We can run the above in Python console：

from option_pricer import EuropeanVanillaPricer
pricer = EuropeanVanillaPricer()
pricer.getMCPrice()
2.1620364099067015
pricer.getBSPrice()
2.1736062697012564

As we can see in the above Monte Carlo simulation, we rely on drawing random numbers, $\epsilon$ from a Standard Normal distribution.
Alternatively, we can use random numbers from a Uniform distribution, i.e. equal probability of each random number.

To do this, we combine step 3, 4 and 5, the current option price is obtained by integrating the terminal payoff under the risk neutral measure:

$e^{-rT}\int_{-\infty}^{\infty}g(S_{0}e^{(\mu-\frac{\sigma ^{2}}{2})T + \sigma \sqrt{T} \epsilon}) \frac{1}{\sqrt{2\pi}}e^{-\frac{\epsilon ^{2}}{2}}d\epsilon$

$=\int_{-\infty}^{\infty}h(\epsilon)\frac{1}{\sqrt{2\pi}}e^{-\frac{\epsilon ^{2}}{2}}d\epsilon$

$=\int_{0}^{1}h[\phi ^{-1}(x)]dx$

$=\int_{0}^{1}f(x)dx$

In the first line, function $g$ is just the payoff condition at expiry. As we are integrating with regard to $\epsilon$ , which follows Standard Normal distribution, the last term is the probability density function.

In the second line, we just use a new function h of epsilon to make the expression more compact.

In the third line, we do inverse transformation to integrate with regard to the cumulative probability, $x$ .

So now it becomes an integral of function $f$ over the UNIFORM distribution with range [0, 1].

Now our simulation task becomes taking random number $x$ from the Uniform distribution [0, 1], and then calculate integral of function $f$ using Monte Carlo.

To be more specific, our task has been changed from calculating

$\frac{1}{N}\sum_{i=1}^{N}h(\epsilon _{i}), where \, \epsilon _{i} \sim \mathcal{N}(0, 1)$

To evaluating

$=\int_{0}^{1}f(x)dx$

When we use Monte Carlo to estimate function integral, we may run into problem of random number clustering, which essentially leads to Convergence rate of $O(N^{0.5})$ .

To conquer this issue, instead of using pseudo-random numbers, we can use a deterministic sequence, whose numbers are more equally spaced. And this is exactly what Quasi Monte Carlo (QMC) does. More details on the low-discrepancy sequences can be found in this post.

European Vanilla Option Pricing – Monte Carlo Methods

Leave a Reply Cancel reply