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

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


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

• Thanks a lot, I will study your answer. But have first added the Monte Carlo code and made a few small clarifications in 'F' and 'interest' in the questions. 'interest' is the interest charged on the Financing Level F. And F = K. Aug 2 '20 at 6:05
• To clarify the 'ratio', this is just used to divide the spot price of expensive stocks. It is just makes the spot price of the options a bit more accessible for smaller investors. Aug 2 '20 at 6:23
• I have checked your answer, thanks a lot for this. Still, my question remains. Since (mu, sigma) of the stock are (0.2469, 0.3781), I would expect the future stock price in 12 months to be much higher than the spot price of S = 425. Hence, I would also expect the option price in 12 months to be much higher than the current spot price. But I am not finding that using Hull's formula. So there must be an error somewhere. Aug 2 '20 at 6:24
• I've added to my answer to show how you can price the option in your MC script and return an option price that matches the Hull formula Aug 2 '20 at 6:34
• Thanks a lot, this was really helpful! Aug 2 '20 at 7:28