diff --git a/Python/examples/european-option-adjoint.py b/Python/examples/european-option-adjoint.py new file mode 100644 index 00000000..c0f1d0ca --- /dev/null +++ b/Python/examples/european-option-adjoint.py @@ -0,0 +1,235 @@ +# --- +# jupyter: +# jupytext: +# formats: py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.4.2 +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# %% [markdown] +# # European options with sensitivities +# +# Copyright (©) 2004, 2005, 2006, 2007 StatPro Italia srl +# Copyright (©) 2024 Xcelerit Computing Limited. +# +# This file is part of QuantLibAAD, a Python wrapper for QuantLib enabled +# for risk computation using automatic differentiation. It uses XAD, +# a fast and comprehensive C++ library for automatic differentiation. +# +# QuantLibAAD and XAD are free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# QuantLibAAD is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . +# +# QuantLib is free software: you can redistribute it and/or modify it under the +# terms of the QuantLib license. You should have received a copy of the +# license along with this program; if not, please email +# . The license is also available online at +# . +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the license for more details. + +# %% [markdown] +# This is a port of the `AdjointEuropeanEquityOption` C++ example. It prices +# a portfolio of European put options (several underlyings, each with a range +# of strikes) over a zero curve and a black variance curve, and computes the +# sensitivities of the total portfolio value to all market inputs (the zero +# rates, the volatilities, the strikes, the underlyings, and the dividend +# yield) in a single adjoint pass. The results are cross-checked against +# bump-and-reprice finite differences. + +# %% +import QuantLibAAD as ql +from xad.adj_1st import Tape +from xad import value +import pandas as pd + +# %% [markdown] +# ### Global parameters + +# %% +calendar = ql.TARGET() +todaysDate = ql.Date(15, ql.May, 1998) +settlementDate = ql.Date(17, ql.May, 1998) +ql.Settings.instance().evaluationDate = todaysDate +dayCounter = ql.Actual365Fixed() +maturity = ql.Date(17, ql.May, 1999) + +# setup and activate the tape +tape = Tape() +tape.activate() + +# %% +interactive = "get_ipython" in globals() + +# %% [markdown] +# ### Market data +# +# The zero curve and the variance curve share the same node dates: the +# settlement date plus the offsets below. The first rate is repeated at the +# first offset date, as in the C++ example. + +# %% +offsets = [13, 41, 75, 165, 256, 345, 524, 703] +curveRates = [0.0357, 0.0349, 0.0341, 0.0355, 0.0359, 0.0368, 0.0386, 0.0401] +curveVols = [0.20, 0.18, 0.178, 0.183, 0.192, 0.203, 0.215, 0.208] + +dates = [settlementDate] + [settlementDate + n for n in offsets] + +rates_in = [ql.Real(curveRates[0])] + [ql.Real(r) for r in curveRates] +vols_in = [ql.Real(v) for v in curveVols] +strikes_in = [ql.Real(float(s)) for s in range(10, 80)] +underlyings_in = [ql.Real(float(s)) for s in range(15, 61, 5)] +dividendYield_in = ql.Real(0.01) + +tape.registerInputs(rates_in) +tape.registerInputs(vols_in) +tape.registerInputs(strikes_in) +tape.registerInputs(underlyings_in) +tape.registerInput(dividendYield_in) + + +# %% [markdown] +# ### Pricing +# +# The whole portfolio is priced inside the recording, so that the curves and +# instruments are built from the registered inputs. + + +# %% +def pricePortfolio(rates, vols, strikes, underlyings, dividendYield): + termStructure = ql.ZeroCurve(dates, rates, dayCounter) + dividendTS = ql.FlatForward(settlementDate, dividendYield, dayCounter) + # the variance curve drops the first (settlement) date + volTS = ql.BlackVarianceCurve(settlementDate, dates[1:], vols, dayCounter) + + exercise = ql.EuropeanExercise(maturity) + value = ql.Real(0.0) + for underlying in underlyings: + process = ql.BlackScholesMertonProcess( + ql.QuoteHandle(ql.SimpleQuote(underlying)), + ql.YieldTermStructureHandle(dividendTS), + ql.YieldTermStructureHandle(termStructure), + ql.BlackVolTermStructureHandle(volTS), + ) + engine = ql.AnalyticEuropeanEngine(process) + for strike in strikes: + option = ql.VanillaOption( + ql.PlainVanillaPayoff(ql.Option.Put, strike), exercise + ) + option.setPricingEngine(engine) + value = value + option.NPV() + return value + + +# %% +tape.newRecording() +portfolioValue = pricePortfolio( + rates_in, vols_in, strikes_in, underlyings_in, dividendYield_in +) + +tape.registerOutput(portfolioValue) +portfolioValue.derivative = 1.0 +tape.computeAdjoints() + +print(f"Portfolio value: {portfolioValue:.5f}") + +# %% [markdown] +# ### Sensitivities + +# %% +results = [] +for d, rate in zip(dates, rates_in): + results.append(("rho", f"rate {d}", rate.derivative)) +for d, vol in zip(dates[1:], vols_in): + results.append(("vega", f"vol {d}", vol.derivative)) +for strike in strikes_in: + results.append(("strike sensitivity", f"strike {value(strike):.0f}", strike.derivative)) +for underlying in underlyings_in: + results.append(("delta", f"underlying {value(underlying):.0f}", underlying.derivative)) +results.append(("dividend rho", "dividend yield", dividendYield_in.derivative)) + +# %% +df = pd.DataFrame(results, columns=["Sensitivity", "Input", "Value"]) +df.style.hide(axis="index") + +# %% +if not interactive: + print(df.to_string(index=False)) + +# %% [markdown] +# ### Cross-check by bumping +# +# Each input is bumped in turn and the portfolio repriced without the tape; +# the finite differences should match the adjoints above. Only the largest +# sensitivities are shown. + +# %% +eps = 1e-5 + + +def bumpedValue(which, index): + rates = [value(r) for r in rates_in] + vols = [value(v) for v in vols_in] + strikes = [value(s) for s in strikes_in] + underlyings = [value(u) for u in underlyings_in] + dividendYield = value(dividendYield_in) + if which == "rate": + rates[index] += eps + elif which == "vol": + vols[index] += eps + elif which == "strike": + strikes[index] += eps + elif which == "underlying": + underlyings[index] += eps + else: + dividendYield += eps + return pricePortfolio( + [ql.Real(r) for r in rates], + [ql.Real(v) for v in vols], + [ql.Real(s) for s in strikes], + [ql.Real(u) for u in underlyings], + ql.Real(dividendYield), + ) + + +# %% +checks = [] +named = ( + [("rate", i, f"rate {dates[i]}", r.derivative) for i, r in enumerate(rates_in)] + + [("vol", i, f"vol {dates[i + 1]}", v.derivative) for i, v in enumerate(vols_in)] + + [("strike", i, f"strike {value(s):.0f}", s.derivative) for i, s in enumerate(strikes_in)] + + [("underlying", i, f"underlying {value(u):.0f}", u.derivative) for i, u in enumerate(underlyings_in)] + + [("dividend", 0, "dividend yield", dividendYield_in.derivative)] +) +for which, index, name, adjoint in sorted(named, key=lambda t: -abs(t[3]))[:8]: + fd = (bumpedValue(which, index) - value(portfolioValue)) / eps + checks.append((name, adjoint, fd)) + +# %% +checkDf = pd.DataFrame(checks, columns=["Input", "Adjoint", "Bump and reprice"]) +checkDf.style.hide(axis="index") + +# %% +if not interactive: + print(checkDf.to_string(index=False)) + +# %% +tape.deactivate()