Pricing a double barrier option using Monte Carlo (C++ & Python code included)

Question

I'm trying to price an option with upper and lower barriers using MC where the payoff is $B_u$ when $S_t > B_u$, $B_l$ when $S_t < B_l$ and $S_t$ when $B_l < S_t < B_u$.

I have written code in both Python and C++, each results in the same price but it doesn’t seem intuitively correct. For the parameters below, price = 109.991. If anyone has any pointers to where the error might be/an analytical solution I'd really appreciate it!

C++:

#include <iostream>
#include <random>
#include <math.h>

// Initialize variables
double s0 = 100;          // Price
double vol = 0.4;         // Volatility
double r = 0.01;          // Interest Rate
double t_ = 255;          // Year
int days = 2;             // Days
int N = pow(10,6);        // Simulations
double b_u = 110;         // Upper Barrier (Rebate)
double b_l = 90;          // Lower Barrier (Rebate)

using namespace std;

std::default_random_engine generator;

double asset_price(double p,double vol,int periods)
{
    double mean = 0.0;
    double stdv = 1.0;

    std::normal_distribution<double> distribution(mean,stdv);

    for(int i=0; i < periods; i++)
    {
        double w = distribution(generator);
        p += s0 * exp((r - 0.5 * pow(vol,2)) * days + vol * sqrt(days) * w);
    }
    return p;
}

int main()
{
    // Monte Carlo Payoffs
    double avg = 0.0;

    for(int j=0; j < N; j++)
    {
        double temp = asset_price(s0,vol,days);
        if(temp > b_u)
        {
            double payoff = b_u;
            payoff = payoff * exp(-r/t_ * days);
            avg += payoff;
        }
        else if(temp < b_l)
        {
            double payoff = b_l;
            payoff = payoff * exp(-r/t_ * days);
            avg += payoff;
        }
        else
        {
            double payoff = temp;
            payoff = payoff * exp(-r/t_ * days);
            avg += payoff;
        }
    }

    // Average Payoff Vector
    double price = avg/(double)N;

    // Results
    cout << "MONTE CARLO BARRIER OPTION PRICING" << endl;
    cout << "----------------------------------" << endl;
    cout << "Option price: " << price << endl;
    cout << "Price at t=0: " << s0 << endl;
    cout << "Volatility: " << vol*100 << "%" << endl;
    cout << "Number of simulations: " << N << endl;

    return 0;
}

Python:

import numpy as np
from math import *


def asset_price(p, v, periods):
    w = np.random.normal(0, 1, size=periods)
    for i in range(periods):
        p += s0 * exp((r - 0.5 * v**2) * days + v * sqrt(days) * w[i])
    return p


# Parameters
s0 = 100  # Price
v = 0.4  # Vol
t_ = 255  # Year
r = 0.01  # Interest Rate
days = 2  # Days until option expiration
N = 100000  # Simulations
avg = 0

# Simulation loop
for i in range(N):
    B_U = 110  # Upper barrier
    B_L = 90  # Lower barrier
    temp = asset_price(s0, v, days)
    if temp > B_U:
        payoff = B_U
        payoff = payoff * np.exp(-r / t_ * days)
        avg += payoff
    elif temp < B_L:
        payoff = B_L
        payoff = payoff * np.exp(-r / t_ * days)
        avg += payoff
    else:
        payoff = temp
        payoff = payoff * np.exp(-r / t_ * days)
        avg += payoff

# Average payoffs vector
price = avg / float(N)

# Results
print "MONTE CARLO BARRIER OPTION PRICING"
print "----------------------------------"
print "Option price: ", price
print "Price at t=0: ", s0
print "Volatility: ", v * 100, "%"

Answer 1

Here are at least three mistakes in your code:

1. p += s0 * exp(...) should be p *= exp(...).

2. Your volatility and rates are per annum, so divide the days by 365 (or 255) in your function asset_price.

3. In asset_price you multiply by days inside the loop. However, the loop is already iterating over the days - so you don't take two steps of one day but two steps of two days in your example.

Some more suggestions/remarks:

1. When you implement a pricer, it is always helpful to check edge cases/limiting behavior. In your case for example you could let B_L = 0, B_H = 1000 (high) and the price of your contract should just be the spot. This already fails in your original code.

2. Using the global variable days or r in your function asset_price is bad style. The function arguments should represent its interface. You however mix passing s0 into the argument p but just access the global variable r from within the function.

3. Note that you are pricing a European barrier option where the barrier is only active at maturity. Not sure if that was your intention.

4. You simulate the asset price at every day from now to maturity. As you are only interested in the price at maturity (see previous point), you could just simulate the value in one larger step. I.e. the loop in asset_price is redundant.

---

Here is how the corrected asset_price could look like:

def asset_price(spot, vola, rate, period_count, delta_t):
    w = np.random.normal(0, 1, size=period_count)
    for i in range(period_count):
        spot *= exp((rate - 0.5 * vola**2) * delta_t + vola * sqrt(delta_t) * w[i])
    return spot

You would pass 1.0 / 255.0 as delta_t such that this is the time step in years. Again, the loop in this function is actually completely redundant (see remark 4).