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

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, "%"

• Since the boundary is only active at expiry you should have a fairly simple analytical solution. – Freelunch Feb 22 '18 at 10:12
• If it can be of any help, I wrote my dissertation on pricing barrier options a couple of years ago. The dissertation and the code (matlab) that goes with is publicly available here: github.com/torbonde/dissertation. Note for instance, that in paragraph 1.2.1 I give analytical expressions for barrier options in the one-dimensional Black-Scholes case. I also consider different ways of pricing barrier options, and from these I would recommend using the Sequential Monte Carlo approach. – torbonde Feb 22 '18 at 11:31

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).

• Why does volatility have to be scaled by the days and not the square root of days? – CodesInChaos Feb 21 '18 at 20:22
• That's not what I am saying. Scaling by square root of time is generally correct. As your loop is over the days already, you need v * sqrt(1 / 255) * w[i]. There should be no days in this line. – LocalVolatility Feb 21 '18 at 20:34
• Great, thank you very much for your suggestions @LocalVolatility. – AlexAbrahams Feb 22 '18 at 14:44