Why does the closed formula result for a Barrier option price deviate so strongly from the Monte Carlo approximation?

Question

I am trying to price a down-and-out, leveraged Barrier option using the closed form formula of Hull (2015). When the price of the underlying asset falls and hits a certain barrier (H), the contract becomes worthless. The issuer of these "Turbo warrants" indicates that the option spot price is calculated as follows: $$P = \frac{S - F}{ratio} $$

(Note: the ratio is used in case the price of the underlying asset is high, like in the case of Amazon stock which is around $3,000. The ratio is often 10 or 100.) The structure of the leveraged Barrier is as follows:

However, the result I get deviates strongly from the result I get when I do an approximation using a Monte Carlo simulation of the stock price. Am I implementing the pricing formula wrongly?

Hull's formula looks like this:

My code looks as follows:

import numpy as np
import scipy.stats as si

# Basic input
S = 425.04 # Spot price of the underlying
r = 0.0054 # Risk-free rate 0.54%
T = 252/252 # life of the option: number of trading days until maturity/ total number of trading days
sigma = 0.3781 # Stock price volatility log returns
mu = 0.2469 # Average stock returns
q = 0 # dividend payout
P = 17.25 # This is the current option spot price 

interest = 0.02 # Annual interest charged on the Financing level. 
# This is charged every day by raising F and raising the Barrier (H)
K = 245.84 * (1 + r)**T # Hull uses the variable K, which is the Financing
# level F in the bar chart I added. I use the Future Value of the Financing level 
# as the exercise price, because if the option hits the Barrier 
# you receive the spot price -/- financing level as payout
H = 270.42 # This is the Barrier level
ratio = 10

# Calculate components of the Hull formula for Barriers
Lambda = (r - q + ((sigma**2)/2))/sigma**2
x1 = (np.log(S/H) / (sigma * np.sqrt(T))) + (Lambda * sigma * np.sqrt(T))
y1 = (np.log(H/S) / (sigma * np.sqrt(T))) + (Lambda * sigma * np.sqrt(T))
x2 = x1 - (sigma * np.sqrt(T))
y2 = y1 - (sigma * np.sqrt(T))

N1 = si.norm.cdf(x1, 0.0, 1.0)
N2 = si.norm.cdf(x2, 0.0, 1.0)
N3 = si.norm.cdf(y1, 0.0, 1.0)
N4 = si.norm.cdf(y2, 0.0, 1.0)

W = (S * N1 * np.exp(-q*T))
X = (K * (np.exp(-r*T)) * N2)
Y = (S * (np.exp(-q*T))) * ((H/S)**(2*Lambda)) * N3
Z = (K * (np.exp(-r*T)) * ((H/S)**((2*Lambda)-2)) * N4)

# Implement the Hull formula for pricing of call Barriers at time T
c =  (W - X - Y + Z) / (ratio)
print("The calculated future price of the option is:", c)
print("The difference between the calculated future and the actual current price is:", c - P)

Based on the above Hull implementation, I arrive at a future price for the option of $17.25. This is only slightly above the current option spot price of $17.92.

However, using on a Monte Carlo simulation of the stock price, I arrive at a future stock price of $544.84. Based on that I would expect an option price of around $29 (based on (S - F) / ratio).

Am I implementing Hull's formula incorrectly or are there other explanations for the large deviation between the closed formula and the Monte Carlo approximation?

EDIT

Here is the Monte Carlo code I am using:

import matplotlib.pyplot as plt
from math import log, e

plt.figure(figsize=(15,10))
num_reps = 1000
delta_t = 1/365
steps = int(T/delta_t)


# get the paths and store them in a list of lists
closing_prices = []
def monti_paths(num_reps, S, steps, mu, sigma, delta_t):
    paths = []
    for j in range(num_reps):
        price_path = [S]
        st = S
        drift = (mu - 0.5 * sigma**2) * delta_t
        sigmasqrtt = sigma * np.sqrt(delta_t)
        for i in range(int(steps)):
            diffusion = sigmasqrtt * np.random.normal(0, 1)
            st *= e**(drift + diffusion)
            price_path.append(st)
        paths.append(price_path)
        closing_prices.append(price_path[-1])
        plt.plot(price_path)
    plt.ylabel('stock price',fontsize=15)
    plt.xlabel('steps',fontsize=15)
    plt.axhline(y = S, color = 'r', linestyle = '-') # print latest price TW
    plt.show()
    return paths

a = monti_paths(num_reps, S, steps, mu, sigma, delta_t)
mean_end_price = np.mean(closing_prices)
print("Expected price in 12 months: $", str(round(mean_end_price,2)))

Answer 1

ADDITIONAL:

There are two problems with your MC code:

1. You're simulating with a drift of $\mu$, which is the real-world dynamics - you will only be able to price the option properly in the risk neutral measure, using a drift of $r$ instead
2. You are not pricing the option, but the underlying stock, and you calculate the final price fo that instead of the option payoff. I've added a snippet below that prices the option along each path, and you can see we are now getting roughly a price of 170 for the option payoff

import matplotlib.pyplot as plt
from math import log, e

plt.figure(figsize=(15,10))
num_reps = 1000
delta_t = 1/365
steps = int(T/delta_t)


# get the paths and store them in a list of lists
closing_prices = []
prices_on_path = []
def price_option_on_path(path, strike, barrier):
    if min(path) < barrier:
        return 0
    else:
        return max(0, path[-1] - strike)

def monti_paths(num_reps, S, steps, mu, sigma, delta_t, strike, barrier):
    paths = []

    for j in range(num_reps):
        price_path = [S]
        st = S
        drift = (mu - 0.5 * sigma**2) * delta_t
        sigmasqrtt = sigma * np.sqrt(delta_t)
        for i in range(int(steps)):
            diffusion = sigmasqrtt * np.random.normal(0, 1)
            st *= e**(drift + diffusion)
            price_path.append(st)
        paths.append(price_path)
        closing_prices.append(price_path[-1])
        prices_on_path.append(price_option_on_path(price_path, strike, barrier))
        plt.plot(price_path)

    plt.ylabel('stock price',fontsize=15)
    plt.xlabel('steps',fontsize=15)
    plt.axhline(y = S, color = 'r', linestyle = '-') # print latest price TW
    plt.show()
    return paths

a = monti_paths(num_reps, S, steps, r, sigma, delta_t, K, H)
mean_end_price = np.mean(closing_prices)
discounted_payoff = np.mean(prices_on_path) / (1+r)**T
print("Expected price in 12 months: $", str(round(mean_end_price,2)))
print("Option payoff in 12 months $", str(round(discounted_payoff,2)))

ORIGINAL:

I'm not entirely sure about the leverage or about how the ratio term fits in above, but for the most part you are trying to price a Down-and-Out Call, so I'll treat the problem just like that. Also, it would be very helpful to see the MC code you are using, as it's hard to know what you are comparing to without it.

First, a couple of general points:

• For this barrier option, your spot $S=425$, your barrier $H=270$ and your strike $K=247$
• As rates are low, the forward at 1 year is close to spot; and much, much higher than the barrier/strike, so your call is far in-the-money
• Because $H>K$, the option will knock out before it hits the strike, so our payoff will just be $(S_t - K)$ unless we knock out. As we are far above the barrier this is unlikely, so raw price $P \simeq S - K = 177$

I've re-implemented the pricing problem using QuantLib (code below - if you haven't got it, you just need to pip install QuantLib-Python to run the snippets) and using these parameters and an analytic pricer, I get a price of 172.49.

I switch to the Monte Carlo pricer and get a price of about 172 (depending on rng, number of paths etc.) which seems close to the analytic price and the price you stated above - so I conjecture that you have done something wrong in your MC simulation

import QuantLib as ql

# Setting up our universe with a single asset and rates curves
spot = 425.04
vol = 0.3781
rate, divi_rate = 0.0054, 0.0
today = ql.Date(1, 7, 2020)

day_count = ql.Actual365Fixed()
calendar = ql.NullCalendar()

volatility = ql.BlackConstantVol(today, calendar, vol, day_count)
riskFreeCurve = ql.FlatForward(today, rate, day_count)
diviCurve = ql.FlatForward(today, divi_rate, day_count)

flat_ts = ql.YieldTermStructureHandle(riskFreeCurve)
dividend_ts = ql.YieldTermStructureHandle(riskFreeCurve)
flat_vol = ql.BlackVolTermStructureHandle(volatility)

# Creating and pricing the barrier option with analytic engine
expiry_date = ql.Date(1, 7, 2021)
strike = 247.16
barrier_level = 270.42

barrier = ql.BarrierOption(ql.Barrier.DownOut, barrier_level, 0.0, ql.PlainVanillaPayoff(ql.Option.Call, strike), ql.EuropeanExercise(expiry_date))

process = ql.BlackScholesProcess(ql.QuoteHandle(ql.SimpleQuote(spot)), flat_ts, flat_vol)
engine = ql.AnalyticBarrierEngine(process)
barrier.setPricingEngine(engine)

print(barrier.NPV())
# 172.49171908548465

And below with a MC engine:

num_paths = 100000
time_length, steps = 8, 2
rng_type = "pseudorandom" # could use "lowdiscrepancy"

rng = ql.GaussianRandomSequenceGenerator(ql.UniformRandomSequenceGenerator(steps, ql.UniformRandomGenerator()))
seq = ql.GaussianPathGenerator(process, time_length, steps, rng, False)

mc_engine = ql.MCBarrierEngine(process, rng_type, steps, requiredSamples=num_paths)
barrier.setPricingEngine(mc_engine)

print(barrier.NPV())
# 172.42679538243036