def test_european_vanilla_option_usage(self):

        european_exercise = EuropeanExercise(self.maturity)
        european_option = VanillaOption(self.payoff, european_exercise)

        analytic_european_engine = AnalyticEuropeanEngine(self.black_scholes_merton_process)

        european_option.set_pricing_engine(analytic_european_engine)

        self.assertAlmostEqual(3.844308, european_option.net_present_value, 6)

    def test_american_vanilla_option(self):

        american_exercise = AmericanExercise(self.maturity)
        american_option = VanillaOption(self.payoff, american_exercise)

        engine = BaroneAdesiWhaleyApproximationEngine(self.black_scholes_merton_process)

        american_option.set_pricing_engine(engine)

        self.assertAlmostEqual(4.459628, american_option.net_present_value, 6)

    def test_analytic_versus_black(self):

        settlement_date = today()
        self.settings.evaluation_date = settlement_date

        daycounter = ActualActual()

        exercise_date = settlement_date + 6 * Months

        payoff = PlainVanillaPayoff(Put, 30)

        exercise = EuropeanExercise(exercise_date)

        risk_free_ts = flat_rate(0.1, daycounter)
        dividend_ts = flat_rate(0.04, daycounter)

        s0 = SimpleQuote(32.0)

        v0    = 0.05
        kappa = 5.0
        theta = 0.05
        sigma = 1.0e-4
        rho   = 0.0

        process = HestonProcess(
            risk_free_ts, dividend_ts, s0, v0, kappa, theta, sigma, rho
        )

        option = VanillaOption(payoff, exercise)

        engine = AnalyticHestonEngine(HestonModel(process), 144)

        option.set_pricing_engine(engine)

        calculated = option.net_present_value

        year_fraction = daycounter.year_fraction(
            settlement_date, exercise_date
        )

        forward_price = 32 * np.exp((0.1 - 0.04) * year_fraction)
        expected = blackFormula(
            payoff.type, payoff.strike, forward_price,
            np.sqrt(0.05 * year_fraction)
        ) * np.exp(-0.1 * year_fraction)

        tolerance = 2.0e-7

        self.assertAlmostEqual(
            calculated,
            expected,
            delta=tolerance
        )

def heston_pricer(trade_date, options, params, rates, spot):
    """
    Price a list of European options with heston model.

    """

    spot = SimpleQuote(spot)

    risk_free_ts = df_to_zero_curve(rates[nm.INTEREST_RATE], trade_date)
    dividend_ts = df_to_zero_curve(rates[nm.DIVIDEND_YIELD], trade_date)

    process = HestonProcess(risk_free_ts, dividend_ts, spot, **params)

    model = HestonModel(process)
    engine = AnalyticHestonEngine(model, 64)

    settlement_date = pydate_to_qldate(trade_date)

    settings = Settings()
    settings.evaluation_date = settlement_date

    modeled_values = np.zeros(len(options))

    for index, row in options.T.iteritems():

        expiry_date = row[nm.EXPIRY_DATE]
        strike = row[nm.STRIKE]

        option_type = Call if row[nm.OPTION_TYPE] == nm.CALL_OPTION else Put

        payoff = PlainVanillaPayoff(option_type, strike)

        expiry_qldate = pydate_to_qldate(expiry_date)
        exercise = EuropeanExercise(expiry_qldate)

        option = VanillaOption(payoff, exercise)
        option.set_pricing_engine(engine)

        modeled_values[index] = option.net_present_value

    prices = options.filter(items=[nm.EXPIRY_DATE, nm.STRIKE,
                                   nm.OPTION_TYPE, nm.SPOT])
    prices[nm.PRICE] = modeled_values
    prices[nm.TRADE_DATE] = trade_date

    return prices

    def _get_option_npv(self):
        """ Suboptimal getter for the npv.

        FIXME: We currently have to recreate most of the objects because we do not
        expose enough of the QuantLib api.

        """

        # convert datetime object to QlDate
        maturity = QlDate.from_datetime(self.maturity)

        underlyingH = SimpleQuote(self.underlying)

        # bootstrap the yield/dividend/vol curves
        flat_term_structure = FlatForward(
            reference_date = settlement_date,
            forward = self.risk_free_rate,
            daycounter = self.daycounter
        )

        flat_dividend_ts = FlatForward(
            reference_date = settlement_date,
            forward = self.dividend_yield,
            daycounter = self.daycounter
        )

        flat_vol_ts = BlackConstantVol(
            settlement_date, calendar, self.volatility, self.daycounter
        )

        black_scholes_merton_process = BlackScholesMertonProcess(
            underlyingH, flat_dividend_ts, flat_term_structure,flat_vol_ts
        )

        payoff = PlainVanillaPayoff(self.option_type, self.strike)

        european_exercise = EuropeanExercise(maturity)

        european_option = VanillaOption(payoff, european_exercise)

        analytic_european_engine = AnalyticEuropeanEngine(black_scholes_merton_process)

        european_option.set_pricing_engine(analytic_european_engine)

        return european_option.net_present_value

def main():
    # global data
    todays_date = Date(15, May, 1998)
    Settings.instance().evaluation_date = todays_date
    settlement_date = Date(17, May ,1998)

    risk_free_rate = FlatForward(
        reference_date = settlement_date,
        forward        = 0.06,
        daycounter     = Actual365Fixed()
    )

    # option parameters
    exercise = AmericanExercise(
        earliest_exercise_date = settlement_date,
        latest_exercise_date   = Date(17, May, 1999)
    )
    payoff = PlainVanillaPayoff(Put, 40.0)

    # market data
    underlying = SimpleQuote(36.0)
    volatility = BlackConstantVol(todays_date, TARGET(), 0.20, Actual365Fixed())
    dividend_yield = FlatForward(
        reference_date = settlement_date,
        forward        = 0.00,
        daycounter     = Actual365Fixed()
    )

    # report
    header = '%19s' % 'method' + ' |' + \
            ' |'.join(['%17s' % tag for tag in ['value',
                                                'estimated error',
                                                'actual error' ] ])
    print
    print header
    print '-'*len(header)

    refValue = None
    def report(method, x, dx = None):
        e = '%.4f' % abs(x-refValue)
        x = '%.5f' % x
        if dx:
            dx = '%.4f' % dx
        else:
            dx = 'n/a'
        print '%19s' % method + ' |' + \
            ' |'.join(['%17s' % y for y in [x, dx, e] ])

    # good to go

    process = BlackScholesMertonProcess(
        underlying, dividend_yield, risk_free_rate, volatility
    )

    option = VanillaOption(payoff, exercise)

    refValue = 4.48667344
    report('reference value',refValue)

    # method: analytic

    option.set_pricing_engine(BaroneAdesiWhaleyApproximationEngine(process))
    report('Barone-Adesi-Whaley',option.net_present_value)

    # method: finite differences
    time_steps = 801
    grid_points = 800

    option.set_pricing_engine(FDAmericanEngine('CrankNicolson', process,time_steps,grid_points))
    report('finite differences',option.net_present_value)


    print 'This is work in progress.'
    print 'Some pricing engines are not yet interfaced.'

    return

    option.set_pricing_engine(BjerksundStenslandEngine(process))
    report('Bjerksund-Stensland',option.NPV())

    # method: binomial
    timeSteps = 801

    option.setPricingEngine(BinomialVanillaEngine(process,'jr',timeSteps))
    report('binomial (JR)',option.NPV())

    option.setPricingEngine(BinomialVanillaEngine(process,'crr',timeSteps))
    report('binomial (CRR)',option.NPV())

    option.setPricingEngine(BinomialVanillaEngine(process,'eqp',timeSteps))
    report('binomial (EQP)',option.NPV())

    option.setPricingEngine(BinomialVanillaEngine(process,'trigeorgis',timeSteps))
    report('bin. (Trigeorgis)',option.NPV())

    option.setPricingEngine(BinomialVanillaEngine(process,'tian',timeSteps))
    report('binomial (Tian)',option.NPV())

    option.setPricingEngine(BinomialVanillaEngine(process,'lr',timeSteps))
    report('binomial (LR)',option.NPV())

#.........这里部分代码省略.........
    # ++++++++++++++++++++ Description of the option +++++++++++++++++++++++++++++++++++++++
    # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    
    Option_name = "IBM Option"
    maturity = Date(26, Jan, 2013)
    strike = 190
    option_type = Call 

    # Here, as an implementation exemple, we make the test with borth american and european exercise
    europeanExercise = EuropeanExercise(maturity)
    # The emericanExercise need also the settlement date, as his right to exerce the buy or call start at the settlement date!
    #americanExercise = AmericanExercise(settlementDate, maturity)
    americanExercise = AmericanExercise(maturity, settlementDate)
    
    print("**********************************")
    print("Description of the option:		", Option_name)
    print("Date of maturity:     			", maturity)
    print("Type of the option:   			", option_type)
    print("Strike of the option:		    ", strike)



    # ++++++++++++++++++ Description of the discrete dividends
    # INPUT You have to determine the frequece and rates of the discrete dividend. Here is a sollution, but she's not the only one.
    # Last know dividend:
    dividend			= 0.75 #//0.75
    next_dividend_date	= Date(10,Feb,2012)
    # HERE we have make the assumption that the dividend will grow with the quarterly croissance:
    dividendCroissance	= 1.03
    dividendfrequence	= Period(3, Months)
    dividendDates = []
    dividends = []


    d = next_dividend_date
    while d <= maturity:
        dividendDates.append(d)
        dividends.append(dividend)
        d = d + dividendfrequence
        dividend *= dividendCroissance

    print("Discrete dividends				")
    print("Dates				Dividends		")
    for date, div in zip(dividendDates, dividends):
        print(date, "		", div)

    # ++++++++++++++++++ Description of the final payoff 
    payoff = PlainVanillaPayoff(option_type, strike)

    # ++++++++++++++++++ The OPTIONS : (American and European) with their dividends description:
    dividendEuropeanOption = DividendVanillaOption(
        payoff, europeanExercise, dividendDates, dividends
    )
    dividendAmericanOption = DividendVanillaOption(
        payoff, americanExercise, dividendDates, dividends
    )


    # just too test
    europeanOption = VanillaOption(payoff, europeanExercise)
    americanOption =  VanillaOption(payoff, americanExercise)

    # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    # ++++++++++++++++++++ Description of the pricing  +++++++++++++++++++++++++++++++++++++
    # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    
    # For the european options we have a closed analytic formula: The Black Scholes:
    dividendEuropeanEngine = AnalyticDividendEuropeanEngine(bsProcess)

    # For the american option we have make the choice of the finite difference model with the CrankNicolson scheme
    #		this model need to precise the time and space step
    #		More they are greater, more the calul will be precise.
    americanGirdPoints = 600
    americanTimeSteps	= 600
    dividendAmericanEngine = FDDividendAmericanEngine('CrankNicolson', bsProcess,americanTimeSteps, americanGirdPoints)

    # just to test
    europeanEngine = AnalyticEuropeanEngine(bsProcess)
    americanEngine = FDAmericanEngine('CrankNicolson', bsProcess,americanTimeSteps, americanGirdPoints)


    # ++++++++++++++++++++ Valorisation ++++++++++++++++++++++++++++++++++++++++
        
    # Link the pricing Engine to the option
    dividendEuropeanOption.set_pricing_engine(dividendEuropeanEngine)
    dividendAmericanOption.set_pricing_engine(dividendAmericanEngine)
    
    # just	to test
    europeanOption.set_pricing_engine(europeanEngine)
    americanOption.set_pricing_engine(americanEngine)

    # Now we make all the needing calcul	
    # ... and final results
    print("NPV of the European Option with discrete dividends=0:	{:.4f}".format(dividendEuropeanOption.npv))
    print("NPV of the European Option without dividend:		{:.4f}".format(europeanOption.npv))
    print("NPV of the American Option with discrete dividends=0:	{:.4f}".format(dividendAmericanOption.npv))
    print("NPV of the American Option without dividend:		{:.4f}".format(americanOption.npv))
    # just a single test
    print("ZeroRate with a maturity at ", maturity, ": ", \
            riskFreeTS.zero_rate(maturity, dayCounter, Simple))

    def test_smith(self):
        # test against result published in
        # Journal of Computational Finance Vol. 11/1 Fall 2007
        # An almost exact simulation method for the heston model

        settlement_date = today()
        self.settings.evaluation_date = settlement_date

        daycounter = ActualActual()
        timeToMaturity = 4

        exercise_date = settlement_date + timeToMaturity * 365

        c_payoff = PlainVanillaPayoff(Call, 100)

        exercise = EuropeanExercise(exercise_date)

        risk_free_ts = flat_rate(0., daycounter)
        dividend_ts = flat_rate(0., daycounter)

        s0 = SimpleQuote(100.0)

        v0    = 0.0194
        kappa = 1.0407
        theta = 0.0586
        sigma = 0.5196
        rho   = -.6747

        nb_steps_a = 100
        nb_paths = 20000
        seed = 12347

        process = HestonProcess(
            risk_free_ts, dividend_ts, s0, v0, kappa, theta,
            sigma, rho, QUADRATICEXPONENTIAL)

        model = HestonModel(process)

        option = VanillaOption(c_payoff, exercise)

        engine = AnalyticHestonEngine(model, 144)

        option.set_pricing_engine(engine)

        price_fft  = option.net_present_value

        engine = MCVanillaEngine(
              trait='MCEuropeanHestonEngine',
              generator='PseudoRandom',
              process=process,
              doAntitheticVariate=True,
              stepsPerYear=nb_steps_a,
              requiredSamples=nb_paths,
              seed=seed)

        option.set_pricing_engine(engine)
        price_mc = option.net_present_value

        expected = 15.1796
        tolerance = .05

        self.assertAlmostEqual(price_fft, expected, delta=tolerance)
        self.assertAlmostEqual(price_mc, expected, delta=tolerance)

    def test_bates_det_jump(self):
        # this looks like a bug in QL:
        # Bates Det Jump model does not have sigma as parameter, yet
        # changing sigma changes the result!

        settlement_date = today()
        self.settings.evaluation_date = settlement_date

        daycounter = ActualActual()

        exercise_date = settlement_date + 6 * Months

        payoff = PlainVanillaPayoff(Put, 1290)
        exercise = EuropeanExercise(exercise_date)
        option = VanillaOption(payoff, exercise)

        risk_free_ts = flat_rate(0.02, daycounter)
        dividend_ts = flat_rate(0.04, daycounter)

        spot = 1290

        ival = {'delta': 3.6828677022272715e-06,
        'kappa': 19.02581428347027,
        'kappaLambda': 1.1209758060939223,
        'lambda': 0.06524550732595163,
        'nu': -1.8968106563601956,
        'rho': -0.7480898462264719,
        'sigma': 1.0206363887835108,
        'theta': 0.01965384459461113,
        'thetaLambda': 0.028915397380738218,
        'v0': 0.06566800935242285}

        process = BatesProcess(
        risk_free_ts, dividend_ts, SimpleQuote(spot),
        ival['v0'], ival['kappa'],
        ival['theta'], ival['sigma'], ival['rho'],
        ival['lambda'], ival['nu'], ival['delta'])

        model = BatesDetJumpModel(process,
                ival['kappaLambda'], ival['thetaLambda'])

        engine = BatesDetJumpEngine(model, 64)

        option.set_pricing_engine(engine)

        calc_1 = option.net_present_value

        ival['sigma'] = 1.e-6

        process = BatesProcess(
        risk_free_ts, dividend_ts, SimpleQuote(spot),
        ival['v0'], ival['kappa'],
        ival['theta'], ival['sigma'], ival['rho'],
        ival['lambda'], ival['nu'], ival['delta'])

        model = BatesDetJumpModel(process,
                ival['kappaLambda'], ival['thetaLambda'])
        engine = BatesDetJumpEngine(model, 64)

        option.set_pricing_engine(engine)

        calc_2 = option.net_present_value

        if(abs(calc_1-calc_2) > 1.e-5):
            print('calc 1 %f calc 2 %f' % (calc_1, calc_2))
        self.assertNotEqual(calc_1, calc_2)

    def test_bucket_analysis_option(self):
        
        settings = Settings()
        
        calendar = TARGET()
        
        todays_date = Date(15, May, 1998)
        settlement_date = Date(17, May, 1998)
        
        settings.evaluation_date = todays_date

        option_type = Put
        underlying = 40
        strike = 40
        dividend_yield = 0.00
        risk_free_rate = 0.001
        volatility = 0.20
        maturity = Date(17, May, 1999)
        daycounter = Actual365Fixed()
        
        underlyingH = SimpleQuote(underlying)
        
        payoff = PlainVanillaPayoff(option_type, strike)
        
        
        flat_term_structure = FlatForward(
            reference_date = settlement_date,
            forward        = risk_free_rate,
            daycounter     = daycounter
        )
        flat_dividend_ts = FlatForward(
            reference_date = settlement_date,
            forward        = dividend_yield,
            daycounter     = daycounter
        )
        
        flat_vol_ts = BlackConstantVol(
            settlement_date,
            calendar,
            volatility,
            daycounter
        )
        
        black_scholes_merton_process = BlackScholesMertonProcess(
            underlyingH,
            flat_dividend_ts,
            flat_term_structure,
            flat_vol_ts
        )
        
        european_exercise = EuropeanExercise(maturity)
        european_option = VanillaOption(payoff, european_exercise)
        analytic_european_engine = AnalyticEuropeanEngine(
                    black_scholes_merton_process
                )
        
        european_option.set_pricing_engine(analytic_european_engine)
        
        
        ba_eo= bucket_analysis(
                [[underlyingH]], [european_option], [1], 0.50, 1)

        self.assertTrue(2, ba_eo)
        self.assertTrue(type(tuple), ba_eo) 
        self.assertEqual(1, len(ba_eo[0][0]))
        self.assertEqual(-0.4582666150152517, ba_eo[0][0][0])

underlyingH = SimpleQuote(underlying)

# bootstrap the yield/dividend/vol curves
flat_term_structure = FlatForward(reference_date=settlement_date, forward=risk_free_rate, daycounter=daycounter)

flat_dividend_ts = FlatForward(reference_date=settlement_date, forward=dividend_yield, daycounter=daycounter)

flat_vol_ts = BlackConstantVol(settlement_date, calendar, volatility, daycounter)

black_scholes_merton_process = BlackScholesMertonProcess(
    underlyingH, flat_dividend_ts, flat_term_structure, flat_vol_ts
)

payoff = PlainVanillaPayoff(option_type, strike)

european_exercise = EuropeanExercise(maturity)

european_option = VanillaOption(payoff, european_exercise)


method = "Black-Scholes"
analytic_european_engine = AnalyticEuropeanEngine(black_scholes_merton_process)

european_option.set_pricing_engine(analytic_european_engine)

print("today: %s settlement: %s maturity: %s" % (todays_date, settlement_date, maturity))
print("NPV: %f\n" % european_option.net_present_value)


### EOF #######################################################################

    def test_compare_BsmHW_HestonHW(self):
        """
        From Quantlib test suite
        """

        print("Comparing European option pricing for a BSM " +
              "process with one-factor Hull-White model...")

        dc = Actual365Fixed()

        todays_date = today()
        settings = Settings()
        settings.evaluation_date = todays_date
        tol = 1.e-2

        spot = SimpleQuote(100)

        dates = [todays_date + Period(i, Years) for i in range(40)]

        rates = [0.01 + 0.0002 * np.exp(np.sin(i / 4.0)) for i in range(40)]
        divRates = [0.02 + 0.0001 * np.exp(np.sin(i / 5.0)) for i in range(40)]

        s0 = SimpleQuote(100)

        r_ts = ZeroCurve(dates, rates, dc)
        q_ts = ZeroCurve(dates, divRates, dc)

        vol = SimpleQuote(0.25)
        vol_ts = BlackConstantVol(
            todays_date,
            NullCalendar(),
            vol.value, dc)

        bsm_process = BlackScholesMertonProcess(
            spot, q_ts, r_ts, vol_ts)

        variance = vol.value * vol.value
        hestonProcess = HestonProcess(
            risk_free_rate_ts=r_ts,
            dividend_ts=q_ts,
            s0=s0,
            v0=variance,
            kappa=5.0,
            theta=variance,
            sigma=1e-4,
            rho=0.0)

        hestonModel = HestonModel(hestonProcess)

        hullWhiteModel = HullWhite(r_ts, a=0.01, sigma=0.01)

        bsmhwEngine = AnalyticBSMHullWhiteEngine(
            0.0, bsm_process, hullWhiteModel)

        hestonHwEngine = AnalyticHestonHullWhiteEngine(
            hestonModel, hullWhiteModel, 128)

        tol = 1e-5
        strikes = [0.25, 0.5, 0.75, 0.8, 0.9,
                   1.0, 1.1, 1.2, 1.5, 2.0, 4.0]
        maturities = [1, 2, 3, 5, 10, 15, 20, 25, 30]
        types = [Put, Call]

        for option_type in types:
            for strike in strikes:
                for maturity in maturities:
                    maturity_date = todays_date + Period(maturity, Years)

                    exercise = EuropeanExercise(maturity_date)

                    fwd = strike * s0.value * \
                        q_ts.discount(maturity_date) / \
                        r_ts.discount(maturity_date)

                    payoff = PlainVanillaPayoff(option_type, fwd)

                    option = VanillaOption(payoff, exercise)

                    option.set_pricing_engine(bsmhwEngine)
                    calculated = option.npv

                    option.set_pricing_engine(hestonHwEngine)
                    expected = option.npv

                    if ((np.abs(expected - calculated) > calculated * tol) and
                       (np.abs(expected - calculated) > tol)):

                        print("Failed to reproduce npv")
                        print("strike    : %f" % strike)
                        print("maturity  : %d" % maturity)
                        print("type      : %s" % option_type.name) 

                    self.assertAlmostEqual(expected, calculated,
                                            delta=tol)

    def test_compare_bsm_bsmhw_hestonhw(self):

        dc = Actual365Fixed()

        todays_date = today()
        settings = Settings()
        settings.evaluation_date = todays_date
        tol = 1.e-2

        spot = SimpleQuote(100)

        dates = [todays_date + Period(i, Years) for i in range(40)]

        rates = [0.01 + 0.0002 * np.exp(np.sin(i / 4.0)) for i in range(40)]
        divRates = [0.02 + 0.0001 * np.exp(np.sin(i / 5.0)) for i in range(40)]

        s0 = SimpleQuote(100)

        r_ts = ZeroCurve(dates, rates, dc)
        q_ts = ZeroCurve(dates, divRates, dc)

        vol = SimpleQuote(0.25)
        vol_ts = BlackConstantVol(
            todays_date,
            NullCalendar(),
            vol.value, dc)

        bsm_process = BlackScholesMertonProcess(
            spot, q_ts, r_ts, vol_ts)

        payoff = PlainVanillaPayoff(Call, 100)
        exercise = EuropeanExercise(dates[1])

        option = VanillaOption(payoff, exercise)

        analytic_european_engine = AnalyticEuropeanEngine(bsm_process)

        option.set_pricing_engine(analytic_european_engine)
        npv_bsm = option.npv

        variance = vol.value * vol.value
        hestonProcess = HestonProcess(
            risk_free_rate_ts=r_ts,
            dividend_ts=q_ts,
            s0=s0,
            v0=variance,
            kappa=5.0,
            theta=variance,
            sigma=1e-4,
            rho=0.0)

        hestonModel = HestonModel(hestonProcess)

        hullWhiteModel = HullWhite(r_ts, a=0.01, sigma=0.01)

        bsmhwEngine = AnalyticBSMHullWhiteEngine(
            0.0, bsm_process, hullWhiteModel)

        hestonHwEngine = AnalyticHestonHullWhiteEngine(
            hestonModel, hullWhiteModel, 128)

        hestonEngine = AnalyticHestonEngine(hestonModel, 144)
        option.set_pricing_engine(hestonEngine)

        npv_heston = option.npv

        option.set_pricing_engine(bsmhwEngine)
        npv_bsmhw = option.npv

        option.set_pricing_engine(hestonHwEngine)
        npv_hestonhw = option.npv

        print("calculated with BSM: %f" % npv_bsm)
        print("BSM-HW: %f" % npv_bsmhw)
        print("Heston: %f" % npv_heston)
        print("Heston-HW: %f" % npv_hestonhw)

        self.assertAlmostEqual(npv_bsm, npv_bsmhw, delta=tol)
        self.assertAlmostEqual(npv_bsm, npv_hestonhw, delta=tol)

    def test_bsm_hw(self):
        print("Testing European option pricing for a BSM process" +
              " with one-factor Hull-White model...")

        dc = Actual365Fixed()
        todays_date = today()
        maturity_date = todays_date + Period(20, Years)

        settings = Settings()
        settings.evaluation_date = todays_date

        spot = SimpleQuote(100)

        q_ts = flat_rate(todays_date, 0.04, dc)
        r_ts = flat_rate(todays_date, 0.0525, dc)
        vol_ts = BlackConstantVol(todays_date, NullCalendar(), 0.25, dc)

        hullWhiteModel = HullWhite(r_ts, 0.00883, 0.00526)

        bsm_process = BlackScholesMertonProcess(spot, q_ts,
                                                r_ts, vol_ts)

        exercise = EuropeanExercise(maturity_date)

        fwd = spot.value * q_ts.discount(maturity_date) / \
            r_ts.discount(maturity_date)

        payoff = PlainVanillaPayoff(Call, fwd)

        option = VanillaOption(payoff, exercise)

        tol = 1e-8
        corr = [-0.75, -0.25, 0.0, 0.25, 0.75]
        expectedVol = [0.217064577, 0.243995801, 0.256402830,
                       0.268236596, 0.290461343]

        for c, v in zip(corr, expectedVol):
            bsm_hw_engine = AnalyticBSMHullWhiteEngine(c, bsm_process,
                                                       hullWhiteModel)

            option = VanillaOption(payoff, exercise)
            option.set_pricing_engine(bsm_hw_engine)
            npv = option.npv

            compVolTS = BlackConstantVol(todays_date, NullCalendar(),
                                         v, dc)

            bs_process = BlackScholesMertonProcess(spot, q_ts,
                                                   r_ts, compVolTS)
            bsEngine = AnalyticEuropeanEngine(bs_process)

            comp = VanillaOption(payoff, exercise)
            comp.set_pricing_engine(bsEngine)

            impliedVol = comp.implied_volatility(npv, bs_process,
                                                 1e-10, 500,
                                                 min_vol=0.1,
                                                 max_vol=0.4)

            if (abs(impliedVol - v) > tol):
                print("Failed to reproduce implied volatility cor: %f" % c)
                print("calculated: %f" % impliedVol)
                print("expected  : %f" % v)

            if abs((comp.npv - npv) / npv) > tol:
                print("Failed to reproduce NPV")
                print("calculated: %f" % comp.npv)
                print("expected  : %f" % npv)

            self.assertAlmostEqual(impliedVol, v, delta=tol)
            self.assertAlmostEqual(comp.npv / npv, 1, delta=tol)

    def test_bucket_analysis_option(self):

        settings = Settings()

        calendar = TARGET()

        todays_date = Date(15, May, 1998)
        settlement_date = Date(17, May, 1998)

        settings.evaluation_date = todays_date

        option_type = Put
        underlying = 40
        strike = 40
        dividend_yield = 0.00
        risk_free_rate = 0.001
        volatility = SimpleQuote(0.20)
        maturity = Date(17, May, 1999)
        daycounter = Actual365Fixed()

        underlyingH = SimpleQuote(underlying)

        payoff = PlainVanillaPayoff(option_type, strike)


        flat_term_structure = FlatForward(
            reference_date = settlement_date,
            forward        = risk_free_rate,
            daycounter     = daycounter
        )
        flat_dividend_ts = FlatForward(
            reference_date = settlement_date,
            forward        = dividend_yield,
            daycounter     = daycounter
        )

        flat_vol_ts = BlackConstantVol(
            settlement_date,
            calendar,
            volatility,
            daycounter
        )

        black_scholes_merton_process = BlackScholesMertonProcess(
            underlyingH,
            flat_dividend_ts,
            flat_term_structure,
            flat_vol_ts
        )

        european_exercise = EuropeanExercise(maturity)
        european_option = VanillaOption(payoff, european_exercise)
        analytic_european_engine = AnalyticEuropeanEngine(
            black_scholes_merton_process
        )

        european_option.set_pricing_engine(analytic_european_engine)


        delta, gamma = bucket_analysis(
            [underlyingH, volatility], [european_option], shift=1e-4,
            type=Centered)
        self.assertAlmostEqual(delta[0], european_option.delta)
        self.assertAlmostEqual(delta[1], european_option.vega)
        self.assertAlmostEqual(gamma[0], european_option.gamma, 5)