Implied Volatility

iv.py

# response to question asked on the Quantlib-users mailing list: https://sourceforge.net/p/quantlib/mailman/message/37391738/
# based off of http://gouthamanbalaraman.com/blog/american-option-pricing-quantlib-python.html

import QuantLib as ql
import scipy.optimize

maturity_date = ql.Date(15, 1, 2016)
spot_price = 127.62
strike_price = 130
volatility = 0.20
dividend_rate =  0.0163
option_type = ql.Option.Call
risk_free_rate = 0.001
day_count = ql.Actual365Fixed()
calendar = ql.UnitedStates()
calculation_date = ql.Date(8, 5, 2015)
ql.Settings.instance().evaluationDate = calculation_date
payoff = ql.PlainVanillaPayoff(option_type, strike_price)
settlement = calculation_date
am_exercise = ql.AmericanExercise(settlement, maturity_date)
american_option = ql.VanillaOption(payoff, am_exercise)
eu_exercise = ql.EuropeanExercise(maturity_date)
european_option = ql.VanillaOption(payoff, eu_exercise)
spot_handle = ql.QuoteHandle(ql.SimpleQuote(spot_price))
flat_ts = ql.YieldTermStructureHandle(ql.FlatForward(calculation_date, risk_free_rate, day_count))
dividend_yield = ql.YieldTermStructureHandle(ql.FlatForward(calculation_date, dividend_rate, day_count))

def get_bsm_process(volatility):
    flat_vol_ts = ql.BlackVolTermStructureHandle(ql.BlackConstantVol(calculation_date, calendar, volatility, day_count))
    return ql.BlackScholesMertonProcess(spot_handle, dividend_yield, flat_ts, flat_vol_ts)

def get_european_npv(volatility):
    bsm_process = get_bsm_process(volatility)
    analytic_engine = ql.AnalyticEuropeanEngine(bsm_process)
    european_option.setPricingEngine(analytic_engine)
    return european_option.NPV()

def get_american_npv(volatility):
    bsm_process = get_bsm_process(volatility)
    binomial_engine = ql.BinomialVanillaEngine(bsm_process, "crr", 2000)
    american_option.setPricingEngine(binomial_engine)
    return american_option.NPV()

european_npv = get_european_npv(volatility)
american_npv = get_american_npv(volatility)

print(f'european npv: {european_npv}')
print(f'american npv: {american_npv}')

european_iv = scipy.optimize.brentq(lambda iv: get_european_npv(iv) - european_npv, 0.1, 1.0)
american_iv = scipy.optimize.brentq(lambda iv: get_american_npv(iv) - american_npv, 0.1, 1.0)

print(f'european iv:  {european_iv}')
print(f'american iv:  {american_iv}')