QuantLib Python Cookbook PDF

QuantLib Python Cookbook Luigi Ballabio and Goutham Balaraman This book is for sale at http://leanpub.com/quantlibpythoncookbook This version was published on 2022-10-29 This is a Leanpub book. Leanpub empowers authors and publishers with the Lean Publishing process. Lean Publishing is the act of...

QuantLib Python Cookbook Luigi Ballabio and Goutham Balaraman This book is for sale at http://leanpub.com/quantlibpythoncookbook This version was published on 2022-10-29 This is a Leanpub book. Leanpub empowers authors and publishers with the Lean Publishing process. Lean Publishing is the act of publishing an in-progress ebook using lightweight tools and many iterations to get reader feedback, pivot until you have the right book and build traction once you do. © 2014 - 2022 Luigi Ballabio and Goutham Balaraman Tweet This Book! Please help Luigi Ballabio and Goutham Balaraman by spreading the word about this book on Twitter! The suggested hashtag for this book is #quantlib. Find out what other people are saying about the book by clicking on this link to search for this hashtag on Twitter: #quantlib Also By Luigi Ballabio Implementing QuantLib 构建 QuantLib Implementing QuantLib の和訳 Contents A note on Python and C++......................................... ii Code conventions used in this book................................... iii Basics................................................... 1 1. QuantLib basics............................................. 2 2. Instruments and pricing engines.................................. 13 3. Numerical Greeks calculation.................................... 21 4. Market quotes.............................................. 26 5. Term structures and their reference dates............................ 34 6. Pricing over a range of days..................................... 40 7. A note on random numbers and dimensionality........................ 50 Interest-rate curves.................................... 60 8. EONIA curve bootstrapping..................................... 61 9. Euribor curve bootstrapping..................................... 73 10. Constructing a yield curve...................................... 101 11. Dangerous day-count conventions................................. 107 12. Implied term structures........................................ 111 13. Interest-rate sensitivities via zero spread............................. 121 14. A glitch in forward-rate curves................................... 128 CONTENTS Interest-rate models................................... 134 15. Simulating interest rates using Hull White model....................... 135 16. Thoughts on the convergence of Hull-White model Monte Carlo simulations..... 140 17. Short interest rate model calibration................................ 150 18. Par versus indexed coupons..................................... 158 19. Modeling interest rate swaps using QuantLib.......................... 162 20. Caps and floors............................................. 167 Equity models.......................................... 172 21. Valuing European option using the Heston model....................... 173 22. Volatility smile and Heston model calibration.......................... 176 23. Heston model parameter calibration in QuantLib Python & SciPy............. 187 24. Valuing European and American options............................. 199 25. Valuing options on commodity futures using the Black formula.............. 204 26. Defining rho for the Black process................................. 208 27. Using curves with different day-count conventions...................... 213 Bonds................................................... 217 28. Modeling fixed rate bonds...................................... 218 29. Building irregular bonds....................................... 221 30. Valuation of bonds with credit spreads.............................. 230 31. Modeling callable bonds........................................ 234 32. Discount margin calculation..................................... 239 33. Duration of floating-rate bonds................................... 244 34. Treasury futures contracts...................................... 251 35. Mischievous pricing conventions.................................. 258 CONTENTS 36. More mischievous conventions................................... 263 Appendix............................................... 270 Translating QuantLib Python examples to C++............................ 271 CONTENTS i The authors have used good faith effort in preparation of this book, but make no expressed or implied warranty of any kind and disclaim without limitation all responsibility for errors or omissions. No liability is assumed for incidental or consequential damages in connection with or arising out of the use of the information or programs contained herein. Use of the information and instructions in this book is at your own risk. The cover image is in the public domain and available from the New York Public Library¹. The cover font is Open Sans Condensed, released by Steve Matteson² under the Apache License version 2.0³. ¹http://digitalcollections.nypl.org/items/510d47df-335e-a3d9-e040-e00a18064a99 ²https://twitter.com/SteveMatteson1 ³http://www.apache.org/licenses/LICENSE-2.0 A note on Python and C++ The choice of using the QuantLib Python bindings and Jupyter was due to their interactivity, which make it easier to demonstrate features, and to the fact that the platform provides out of the box excellent modules like matplotlib for graphing and pandas for data analysis. This choice might seem to leave C++ users out in the cold. However, it’s easy enough to translate the Python code shown here into the corresponding C++ code. An example of such translation is shown in the appendix. Code conventions used in this book The recipes in this cookbook are written as Jupyter notebooks⁴, and follow their structure: blocks of explanatory text, like the one you’re reading now, are mixed with cells containing Python code (inputs) and the results of executing it (outputs). The code and its output—if any—are marked by In [N] and Out [N], respectively, with N being the index of the cell. You can see an example in the computations below: In : def f(x, y): return x + 2*y In : a = 4 b = 2 f(a, b) Out: 8 By default, Jupyter displays the result of the last instruction as the output of a cell, like it did above; however, print statements can display further results. In : print(a) print(b) print(f(b, a)) Out: 4 2 10 Jupyter also knows a few specific data types, such as Pandas data frames, and displays them in a more readable way: In : import pandas as pd pd.DataFrame({ 'foo': [1,2,3], 'bar': ['a','b','c'] }) Out: ⁴http://jupyter.org/ Code conventions used in this book iv foo bar 0 1 a 1 2 b 2 3 c The index of the cells shows the order of their execution. Jupyter doesn’t constrain it; however, in all of the recipes of this book the cells were executed in sequential order as displayed. All cells are executed in the global Python scope; this means that, as we execute the code in the recipes, all variables, functions and classes defined in a cell are available to the ones that follow. Notebooks can also include plots, as in the following cell: In : %matplotlib inline import numpy as np import utils f, ax = utils.plot(figsize=(10,2)) ax.plot([0, 0.25, 0.5, 0.75, 1.0], np.random.random(5)) Out: [] As you might have noted, the cell above also printed a textual representation of the object returned from the plot, since it’s the result of the last instruction in the cell. To prevent this, cells in the recipes might have a semicolon at the end, as in the next cell. This is just a quirk of the Jupyter display system, and it doesn’t have any particular significance; I’m mentioning it here just so that you dont’t get confused by it. In : f, ax = utils.plot(figsize=(10,2)) ax.plot([0, 0.25, 0.5, 0.75, 1.0], np.random.random(5)); Code conventions used in this book v Finally, the utils module that I imported above is a short module containing convenience functions, mostly related to plots, for the notebooks in this collection. It’s not necessary to understand its implementation to follow the recipes, and therefore we won’t cover it here; but if you’re interested and want to look at it, it’s included in the zip archive that you can download from Leanpub if you purchased the book. Basics 1. QuantLib basics In this chapter we will introduce some of the basic concepts such as Date, Period, Calendar and Schedule. These are QuantLib constructs that are used throughout the library in creation of instruments, models, term structures etc. In : import QuantLib as ql import pandas as pd Date Class The Date object can be created using the constructor as Date(day, month, year). It would be worthwhile to pay attention to the fact that day is the first argument, followed by month and then the year. This is different from the Python datetime object instantiation. In : date = ql.Date(31, 3, 2015) print(date) Out: March 31st, 2015 The fields of the Date object can be accessed using the month(), dayOfMonth() and year() methods. The weekday() method can be used to fetch the day of the week. In : print("%d-%d-%d" %(date.month(), date.dayOfMonth(), date.year())) Out: 3-31-2015 In : date.weekday() == ql.Tuesday Out: True The Date objects can also be used to perform arithmetic operations such as advancing by days, weeks, months etc. Periods such as weeks or months can be denoted using the Period class. Period object constructor signature is Period(num_periods, period_type). The num_periods is an integer and represents the number of periods. The period_type can be Weeks, Months and Years. QuantLib basics 3 In : type(date+1) Out: QuantLib.QuantLib.Date In : print("Add a day : {0}".format(date + 1)) print("Subtract a day : {0}".format(date - 1)) print("Add a week : {0}".format(date + ql.Period(1, ql.Weeks))) print("Add a month : {0}".format(date + ql.Period(1, ql.Months))) print("Add a year : {0}".format(date + ql.Period(1, ql.Years))) Out: Add a day : April 1st, 2015 Subtract a day : March 30th, 2015 Add a week : April 7th, 2015 Add a month : April 30th, 2015 Add a year : March 31st, 2016 One can also do logical operations using the Date object. In : print(date == ql.Date(31, 3, 2015)) print(date > ql.Date(30, 3, 2015)) print(date < ql.Date(1, 4, 2015)) print(date != ql.Date(1, 4, 2015)) Out: True True True True The Date object is used in setting valuation dates, issuance and expiry dates of instruments. The Period object is used in setting tenors, such as that of coupon payments, or in constructing payment schedules. Calendar Class The Date arithmetic above did not take holidays into account. But valuation of different securities would require taking into account the holidays observed in a specific exchange or country. The Calendar class implements this functionality for all the major exchanges. Let us take a look at a few examples here. QuantLib basics 4 In : date = ql.Date(31, 3, 2015) us_calendar = ql.UnitedStates(ql.UnitedStates.GovernmentBond) italy_calendar = ql.Italy() period = ql.Period(60, ql.Days) raw_date = date + period us_date = us_calendar.advance(date, period) italy_date = italy_calendar.advance(date, period) print("Add 60 days: {0}".format(raw_date)) print("Add 60 business days in US: {0}".format(us_date)) print("Add 60 business days in Italy: {0}".format(italy_date)) Out: Add 60 days: May 30th, 2015 Add 60 business days in US: June 24th, 2015 Add 60 business days in Italy: June 26th, 2015 The addHoliday and removeHoliday methods in the calendar can be used to add and remove holidays to the calendar respectively. If a calendar has any missing holidays or has a wrong holiday, then these methods come handy in fixing the errors. The businessDaysBetween method helps find out the number of business days between two dates per a given calendar. Let us use this method on the us_calendar and italy_calendar as a sanity check. In : us_busdays = us_calendar.businessDaysBetween(date, us_date) italy_busdays = italy_calendar.businessDaysBetween(date, italy_date) print("Business days US: {0}".format(us_busdays)) print("Business days Italy: {0}".format(italy_busdays)) Out: Business days US: 60 Business days Italy: 60 In valuation of certain deals, more than one calendar’s holidays are observed. QuantLib has JointCalendar class that allows you to combine the holidays of two or more calendars. Let us take a look at a working example. QuantLib basics 5 In : joint_calendar = ql.JointCalendar(us_calendar, italy_calendar) joint_date = joint_calendar.advance(date, period) joint_busdays = joint_calendar.businessDaysBetween(date, joint_date) print("Add 60 business days in US-Italy: {0}".format(joint_date)) print("Business days US-Italy: {0}".format(joint_busdays)) Out: Add 60 business days in US-Italy: June 29th, 2015 Business days US-Italy: 60 Schedule Class The Schedule object is necessary in creating coupon schedules or call schedules. Schedule object constructors have the following signature: Schedule(const Date& effectiveDate, const Date& terminationDate, const Period& tenor, const Calendar& calendar, BusinessDayConvention convention, BusinessDayConvention terminationDateConvention, DateGeneration::Rule rule, bool endOfMonth, const Date& firstDate = Date(), const Date& nextToLastDate = Date()) and Schedule(const std::vector&, const Calendar& calendar, BusinessDayConvention rollingConvention) In : effective_date = ql.Date(1, 1, 2015) termination_date = ql.Date(1, 1, 2016) tenor = ql.Period(ql.Monthly) calendar = ql.UnitedStates(ql.UnitedStates.GovernmentBond) business_convention = ql.Following termination_business_convention = ql.Following date_generation = ql.DateGeneration.Forward end_of_month = False schedule = ql.Schedule(effective_date, termination_date, tenor, QuantLib basics 6 calendar, business_convention, termination_business_convention, date_generation, end_of_month) pd.DataFrame({'date': list(schedule)}) Out: date 0 January 2nd, 2015 1 February 2nd, 2015 2 March 2nd, 2015 3 April 1st, 2015 4 May 1st, 2015 5 June 1st, 2015 6 July 1st, 2015 7 August 3rd, 2015 8 September 1st, 2015 9 October 1st, 2015 10 November 2nd, 2015 11 December 1st, 2015 12 January 4th, 2016 Here we have generated a Schedule object that will contain dates between effective_date and termination_date with the tenor specifying the Period to be Monthly. The calendar object is used for determining holidays. Here we have chosen the convention to be the day following holidays. That is why we see that holidays are excluded in the list of dates. The Schedule class can handle generation of dates with irregularity in schedule. The two extra parameters firstDate and nextToLastDate parameters along with a combination of forward or backward date generation rule can be used to generate short or long stub payments at the front or back end of the schedule. For example, the following combination of firstDate and backward generation rule creates a short stub in the front on the January 15, 2015. In : # short stub in the front effective_date = ql.Date(1, 1, 2015) termination_date = ql.Date(1, 1, 2016) first_date = ql.Date(15, 1, 2015) schedule = ql.Schedule(effective_date, termination_date, tenor, calendar, business_convention, termination_business_convention, QuantLib basics 7 ql.DateGeneration.Backward, end_of_month, first_date) pd.DataFrame({'date': list(schedule)}) Out: date 0 January 2nd, 2015 1 January 15th, 2015 2 February 2nd, 2015 3 March 2nd, 2015 4 April 1st, 2015 5 May 1st, 2015 6 June 1st, 2015 7 July 1st, 2015 8 August 3rd, 2015 9 September 1st, 2015 10 October 1st, 2015 11 November 2nd, 2015 12 December 1st, 2015 13 January 4th, 2016 Using the nextToLastDate parameter along with the forward date generation rule creates a short stub at the back end of the schedule. In : # short stub at the back effective_date = ql.Date(1, 1, 2015) termination_date = ql.Date(1, 1, 2016) penultimate_date = ql.Date(15, 12, 2015) schedule = ql.Schedule(effective_date, termination_date, tenor, calendar, business_convention, termination_business_convention, ql.DateGeneration.Forward, end_of_month, ql.Date(), penultimate_date) pd.DataFrame({'date': list(schedule)}) Out: QuantLib basics 8 date 0 January 2nd, 2015 1 February 2nd, 2015 2 March 2nd, 2015 3 April 1st, 2015 4 May 1st, 2015 5 June 1st, 2015 6 July 1st, 2015 7 August 3rd, 2015 8 September 1st, 2015 9 October 1st, 2015 10 November 2nd, 2015 11 December 1st, 2015 12 December 15th, 2015 13 January 4th, 2016 Using the backward generation rule along with the firstDate allows us to create a long stub in the front. Below the first two dates are longer in duration than the rest of the dates. In : # long stub in the front first_date = ql.Date(1, 2, 2015) effective_date = ql.Date(15, 12, 2014) termination_date = ql.Date(1, 1, 2016) schedule = ql.Schedule(effective_date, termination_date, tenor, calendar, business_convention, termination_business_convention, ql.DateGeneration.Backward, end_of_month, first_date) pd.DataFrame({'date': list(schedule)}) Out: QuantLib basics 9 date 0 December 15th, 2014 1 February 2nd, 2015 2 March 2nd, 2015 3 April 1st, 2015 4 May 1st, 2015 5 June 1st, 2015 6 July 1st, 2015 7 August 3rd, 2015 8 September 1st, 2015 9 October 1st, 2015 10 November 2nd, 2015 11 December 1st, 2015 12 January 4th, 2016 Similarly the usage of nextToLastDate parameter along with forward date generation rule can be used to generate long stub at the back of the schedule. In : # long stub at the back effective_date = ql.Date(1, 1, 2015) penultimate_date = ql.Date(1, 12, 2015) termination_date = ql.Date(15, 1, 2016) schedule = ql.Schedule(effective_date, termination_date, tenor, calendar, business_convention, termination_business_convention, ql.DateGeneration.Forward, end_of_month, ql.Date(), penultimate_date) pd.DataFrame({'date': list(schedule)}) Out: QuantLib basics 10 date 0 January 2nd, 2015 1 February 2nd, 2015 2 March 2nd, 2015 3 April 1st, 2015 4 May 1st, 2015 5 June 1st, 2015 6 July 1st, 2015 7 August 3rd, 2015 8 September 1st, 2015 9 October 1st, 2015 10 November 2nd, 2015 11 December 1st, 2015 12 January 15th, 2016 Below the Schedule is generated from a list of dates. In : dates = [ql.Date(2,1,2015), ql.Date(2, 2,2015), ql.Date(2,3,2015), ql.Date(1,4,2015), ql.Date(1,5,2015), ql.Date(1,6,2015), ql.Date(1,7,2015), ql.Date(3,8,2015), ql.Date(1,9,2015), ql.Date(1,10,2015), ql.Date(2,11,2015), ql.Date(1,12,2015), ql.Date(4,1,2016)] rolling_convention = ql.Following schedule = ql.Schedule(dates, calendar, rolling_convention) pd.DataFrame({'date': list(schedule)}) Out: date 0 January 2nd, 2015 1 February 2nd, 2015 2 March 2nd, 2015 3 April 1st, 2015 4 May 1st, 2015 5 June 1st, 2015 6 July 1st, 2015 7 August 3rd, 2015 8 September 1st, 2015 9 October 1st, 2015 10 November 2nd, 2015 11 December 1st, 2015 12 January 4th, 2016 QuantLib basics 11 Interest Rate The InterestRate class can be used to store the interest rate with the compounding type, day count and the frequency of compounding. Below we show how to create an interest rate of 5.0% compounded annually, using Actual/Actual day count convention. In : annual_rate = 0.05 day_count = ql.ActualActual(ql.ActualActual.ISDA) compound_type = ql.Compounded frequency = ql.Annual interest_rate = ql.InterestRate(annual_rate, day_count, compound_type, frequency) print(interest_rate) Out: 5.000000 % Actual/Actual (ISDA) Annual compounding Lets say if you invest a dollar at the interest rate described by interest_rate, the compoundFactor method in the InterestRate object gives you how much your investment will be worth after any period. Below we show that the value returned by compound_factor for 2 years agrees with the expected compounding formula. In : t = 2.0 print(interest_rate.compoundFactor(t)) print((1+annual_rate)*(1.0+annual_rate)) Out: 1.1025 1.1025 The discountFactor method returns the reciprocal of the compoundFactor method. The discount factor is useful while calculating the present value of future cashflows. In : print(interest_rate.discountFactor(t)) print(1.0/interest_rate.compoundFactor(t)) Out: 0.9070294784580498 0.9070294784580498 A given interest rate can be converted into other compounding types and compounding frequency using the equivalentRate method. QuantLib basics 12 In : new_frequency = ql.Semiannual new_interest_rate = interest_rate.equivalentRate(compound_type, new_frequency, t) print(new_interest_rate) Out: 4.939015 % Actual/Actual (ISDA) Semiannual compounding The discount factor for the two InterestRate objects, interest_rate and new_interest_rate are the same, as shown below. In : print(interest_rate.discountFactor(t)) print(new_interest_rate.discountFactor(t)) Out: 0.9070294784580498 0.9070294784580495 The impliedRate method in the InterestRate class takes compound factor to return the implied rate. The impliedRate method is a static method in the InterestRate class and can be used without an instance of InterestRate. Internally the equivalentRate method invokes the impliedRate method in its calculations. Conclusion This chapter gave an introduction to the basics of QuantLib. Here we explained the Date, Schedule, Calendar and InterestRate classes. 2. Instruments and pricing engines In this notebook, I’ll show how instruments and their available engines can monitor changes in their input data. Setup To begin, we import the QuantLib module and set up the global evaluation date. In : import QuantLib as ql In : today = ql.Date(7, ql.March, 2014) ql.Settings.instance().evaluationDate = today The instrument As a sample instrument, we’ll take a textbook example: a European option. Building the option requires only the specification of its contract, so its payoff (it’s a call option with strike at 100) and its exercise, three months from today’s date. Market data will be selected and passed later, depending on the calculation methods. In : option = ql.EuropeanOption(ql.PlainVanillaPayoff(ql.Option.Call, 100.0), ql.EuropeanExercise(ql.Date(7, ql.June, 2014))) First pricing method: analytic Black-Scholes formula The different pricing methods are implemented as pricing engines holding the required market data. The first we’ll use is the one encapsulating the analytic Black-Scholes formula. First, we collect the quoted market data. We’ll assume flat risk-free rate and volatility, so they can be expressed by SimpleQuote instances: those model numbers whose value can change and that can notify observers when this happens. The underlying value is at 100, the risk-free value at 1%, and the volatility at 20%. Instruments and pricing engines 14 In : u = ql.SimpleQuote(100.0) r = ql.SimpleQuote(0.01) sigma = ql.SimpleQuote(0.20) In order to build the engine, the market data are encapsulated in a Black-Scholes process object. First we build flat curves for the risk-free rate and the volatility… In : riskFreeCurve = ql.FlatForward(0, ql.TARGET(), ql.QuoteHandle(r), ql.Actual360()) volatility = ql.BlackConstantVol(0, ql.TARGET(), ql.QuoteHandle(sigma), ql.Actual360()) …then we instantiate the process with the underlying value and the curves we just built. The inputs are all stored into handles, so that we could change the quotes and curves used if we wanted. I’ll skip over this for the time being. In : process = ql.BlackScholesProcess(ql.QuoteHandle(u), ql.YieldTermStructureHandle(riskFreeCurve), ql.BlackVolTermStructureHandle(volatility)) Once we have the process, we can finally use it to build the engine… In : engine = ql.AnalyticEuropeanEngine(process) …and once we have the engine, we can set it to the option and evaluate the latter. In : option.setPricingEngine(engine) In : print(option.NPV()) Out: 4.155543462156206 Depending on the instrument and the engine, we can also ask for other results; in this case, we can ask for Greeks. Instruments and pricing engines 15 In : print(option.delta()) print(option.gamma()) print(option.vega()) Out: 0.5302223303784392 0.03934493301271913 20.109632428723106 Market changes As I mentioned, market data are stored in Quote instances and thus can notify the option when any of them changes. We don’t have to do anything explicitly to tell the option to recalculate: once we set a new value to the underlying, we can simply ask the option for its NPV again and we’ll get the updated value. In : u.setValue(105.0) print(option.NPV()) Out: 7.27556357927846 Just for showing off, we can use this to graph the option value depending on the underlying asset value. After a bit of graphic setup (don’t pay attention to the man behind the curtains)… In : %matplotlib inline import numpy as np from IPython.display import display import utils …we can take an array of values from 80 to 120, set the underlying value to each of them, collect the corresponding option values, and plot the results. In : f, ax = utils.plot() xs = np.linspace(80.0, 120.0, 400) ys = [] for x in xs: u.setValue(x) ys.append(option.NPV()) ax.set_title('Option value') utils.highlight_x_axis(ax) ax.plot(xs, ys); Instruments and pricing engines 16 Other market data also affect the value, of course. In : u.setValue(105.0) r.setValue(0.01) sigma.setValue(0.20) In : print(option.NPV()) Out: 7.27556357927846 We can see it when we change the risk-free rate… In : r.setValue(0.03) In : print(option.NPV()) Out: 7.624029148527754 …or the volatility. Instruments and pricing engines 17 In : sigma.setValue(0.25) In : print(option.NPV()) Out: 8.531296969971573 Date changes Just as it does when inputs are modified, the value also changes if we advance the evaluation date. Let’s look first at the value of the option when its underlying is worth 105 and there’s still three months to exercise… In : u.setValue(105.0) r.setValue(0.01) sigma.setValue(0.20) print(option.NPV()) Out: 7.27556357927846 …and then move to a date two months before exercise. In : ql.Settings.instance().evaluationDate = ql.Date(7, ql.April, 2014) Again, we don’t have to do anything explicitly: we just ask the option its value, and as expected it has decreased, as can also be seen by updating the plot. In : print(option.NPV()) Out: 6.560073820974377 In : ys = [] for x in xs: u.setValue(x) ys.append(option.NPV()) ax.plot(xs, ys, '--') display(f) Instruments and pricing engines 18 In the default library configuration, the returned value goes down to 0 when we reach the exercise date. In : ql.Settings.instance().evaluationDate = ql.Date(7, ql.June, 2014) In : print(option.NPV()) Out: 0.0 Other pricing methods The pricing-engine mechanism allows us to use different pricing methods. For comparison, I’ll first set the input data back to what they were previously and output the Black-Scholes price. Instruments and pricing engines 19 In : ql.Settings.instance().evaluationDate = today u.setValue(105.0) r.setValue(0.01) sigma.setValue(0.20) In : print(option.NPV()) Out: 7.27556357927846 Let’s say that we want to use a Heston model to price the option. What we have to do is to instantiate the corresponding class with the desired inputs… In : model = ql.HestonModel( ql.HestonProcess( ql.YieldTermStructureHandle(riskFreeCurve), ql.YieldTermStructureHandle(ql.FlatForward(0, ql.TARGET(), 0.0, ql.Actual360())), ql.QuoteHandle(u), 0.04, 0.1, 0.01, 0.05, -0.75)) …pass it to the corresponding engine, and set the new engine to the option. In : engine = ql.AnalyticHestonEngine(model) option.setPricingEngine(engine) Asking the option for its NPV will now return the value according to the new model. In : print(option.NPV()) Out: 7.295356086978629 Lazy recalculation One last thing. Up to now, we haven’t really seen evidence of notifications going around. After all, the instrument might just have recalculated its value every time, regardless of notifications. What I’m going to show, instead, is that the option doesn’t just recalculate every time anything changes; it also avoids recalculations when nothing has changed. We’ll switch to a Monte Carlo engine, which takes a few seconds to run the required simulation. Instruments and pricing engines 20 In : engine = ql.MCEuropeanEngine(process, "PseudoRandom", timeSteps=20, requiredSamples=250000) option.setPricingEngine(engine) When we ask for the option value, we have to wait for the calculation to finish… In : %time print(option.NPV()) Out: 7.248816340995633 CPU times: user 1.59 s, sys: 0 ns, total: 1.59 s Wall time: 1.58 s …but a second call to the NPV method will be instantaneous when made before anything changes. In this case, the option didn’t calculate its value; it just returned the result that it cached from the previous call. In : %time print(option.NPV()) Out: 7.248816340995633 CPU times: user 1.53 ms, sys: 18 µs, total: 1.55 ms Wall time: 1.65 ms If we change anything (e.g., the underlying value)… In : u.setValue(104.0) …the option is notified of the change, and the next call to NPV will again take a while. In : %time print(option.NPV()) Out: 6.598508180782789 CPU times: user 1.56 s, sys: 4.84 ms, total: 1.57 s Wall time: 1.56 s 3. Numerical Greeks calculation In this notebook, I’ll build on the facilities provided by the Instrument class (that is, its ability to detect changes in its inputs and recalculate accordingly) to show how to calculate numerical Greeks when the engine doesn’t provide them. Setup As usual, we import the QuantLib module and set the evaluation date: In : import QuantLib as ql In : today = ql.Date(8, ql.October, 2014) ql.Settings.instance().evaluationDate = today A somewhat exotic option As an example, we’ll use a knock-in barrier option: In : option = ql.BarrierOption(ql.Barrier.UpIn, 120.0, # barrier 0.0, # rebate ql.PlainVanillaPayoff(ql.Option.Call, 100.0), ql.EuropeanExercise(ql.Date(8, ql.January, 2015))) For the purpose of this example, the market data are the underlying value, the risk-free rate and the volatility. We wrap them in quotes, so that the instrument will be notified of any changes… In : u = ql.SimpleQuote(100.0) r = ql.SimpleQuote(0.01) sigma = ql.SimpleQuote(0.20) …and from the quotes we build the flat curves and the process that the engine requires. As explained in a later notebook, we build the term structures so that they move with the evaluation date; this will be useful further on. Numerical Greeks calculation 22 In : riskFreeCurve = ql.FlatForward(0, ql.TARGET(), ql.QuoteHandle(r), ql.Actual360()) volatility = ql.BlackConstantVol(0, ql.TARGET(), ql.QuoteHandle(sigma), ql.Actual360()) In : process = ql.BlackScholesProcess(ql.QuoteHandle(u), ql.YieldTermStructureHandle(riskFreeCurve), ql.BlackVolTermStructureHandle(volatility)) Finally, we build the engine (the library provides one based on an analytic formula) and set it to the option. In : option.setPricingEngine(ql.AnalyticBarrierEngine(process)) Now we can ask the option for its value… In : print(option.NPV()) Out: 1.3657980739109867 …but we’re not so lucky when it comes to Greeks: In : print(option.delta()) Out: --------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) Cell In , line 1 ----> 1 print(option.delta()) File /usr/local/lib/python3.10/site-packages/QuantLib/QuantLib.py:11925, in OneAssetOp\ tion.delta(self) 11924 def delta(self): > 11925 return _QuantLib.OneAssetOption_delta(self) RuntimeError: delta not provided The engine doesn’t provide the delta, so asking for it raises an error. Numerical calculation What does a quant have to do? We can use numerical differentiation to approximate the Greeks, as shown in the next figure: that is, we can approximate the derivative by calculating the option value for two slightly different values of the underlying and by taking the slope between the resulting points. Numerical Greeks calculation 23 The relevant formulas are: P (u0 + h) − P (u0 − h) P (u0 + h) − 2P (u0 ) + P (u0 − h) ∆= Γ= 2h h2 where P (u) is the price of the option for a given value of the underlying u. Thanks to the framework we set in place, getting the perturbed prices is easy enough. We just have to set the relevant quote to the new value and ask the option for its price again. Thus, we choose a small increment and start. First, we save the current value of the option… In : u0 = u.value() ; h = 0.01 In : P0 = option.NPV() ; print(P0) Out: 1.3657980739109867 …then we increase the underlying value and get the new option value… In : u.setValue(u0+h) P_plus = option.NPV() ; print(P_plus) Out: 1.3688112201958083 …then we do the same after decreasing the underlying value. Numerical Greeks calculation 24 In : u.setValue(u0-h) P_minus = option.NPV() ; print(P_minus) Out: 1.3627900998610207 Finally, we set the underlying value back to its current value. In : u.setValue(u0) Applying the formulas above give us the desired Greeks: In : Delta = (P_plus - P_minus)/(2*h) Gamma = (P_plus - 2*P0 + P_minus)/(h*h) print(Delta) print(Gamma) Out: 0.3010560167393761 0.05172234855521651 The approach is usable for any Greek. We can use the two-sided formula above, or the one-sided formula below if we want to minimize the number of evaluations: ∂P P (x0 + h) − P (x0 ) = ∂x h For instance, here we calculate Rho and Vega: In : r0 = r.value() ; h = 0.0001 r.setValue(r0+h) ; P_plus = option.NPV() r.setValue(r0) Rho = (P_plus - P0)/h ; print(Rho) Out: 6.531038494277386 In : sigma0 = sigma.value() ; h = 0.0001 sigma.setValue(sigma0+h) ; P_plus = option.NPV() sigma.setValue(sigma0) Vega = (P_plus - P0)/h ; print(Vega) Out: 26.52519924198904 The approach for the Theta is a bit different, although it still relies on the fact that the option reacts to the change in the market data. The problem is that we don’t have the time to maturity available as a quote, as was the case for the other quantities. Instead, since we set up the term structures so that they move with the evaluation date, we just have to set it to tomorrow’s date to get the corresponding option value: Numerical Greeks calculation 25 In : ql.Settings.instance().evaluationDate = today+1 P1 = option.NPV() h = 1.0/365 Theta = (P1-P0)/h ; print(Theta) Out: -10.770888399441302 4. Market quotes In this notebook, I’ll show a pitfall to avoid when multiple quotes need to be updated. In : %matplotlib inline import numpy as np import utils In : import QuantLib as ql In : today = ql.Date(17, ql.October, 2016) ql.Settings.instance().evaluationDate = today Setting the stage For illustration purposes, I’ll create a bond curve using the same data and algorithm shown in one of the QuantLib C++ examples; namely, I’ll give to the curve the functional form defined by the Nelson-Siegel model and I’ll fit it to a number of bond. Here are the maturities in years and the coupons of the bonds I’ll use: In : data = [ (2, 0.02), (4, 0.0225), (6, 0.025), (8, 0.0275), (10, 0.03), (12, 0.0325), (14, 0.035), (16, 0.0375), (18, 0.04), (20, 0.0425), (22, 0.045), (24, 0.0475), (26, 0.05), (28, 0.0525), (30, 0.055)] For simplicity, I’ll use the same start date, frequency and conventions for all the bonds; this doesn’t affect the point I’m going to make in the rest of the notebook. I’ll also assume that all bonds currently price at 100. I’ll skip over the details of building the curve now; the one thing you’ll need to remember is that it depends on the quotes modeling the bond prices. Market quotes 27 In : calendar = ql.TARGET() settlement = calendar.advance(today, 3, ql.Days) quotes = [] helpers = [] for length, coupon in data: maturity = calendar.advance(settlement, length, ql.Years) schedule = ql.Schedule( settlement, maturity, ql.Period(ql.Annual), calendar, ql.ModifiedFollowing, ql.ModifiedFollowing, ql.DateGeneration.Backward, False) quote = ql.SimpleQuote(100.0) quotes.append(quote) helpers.append( ql.FixedRateBondHelper(ql.QuoteHandle(quote), 3, 100.0, schedule, [coupon], ql.SimpleDayCounter(), ql.ModifiedFollowing)) curve = ql.FittedBondDiscountCurve(0, calendar, helpers, ql.SimpleDayCounter(), ql.NelsonSiegelFitting()) Here is a visualization of the curve as discount factors versus time in years: In : sample_times = np.linspace(0.0, 30.0, 301) sample_discounts = [ curve.discount(t) for t in sample_times ] f, ax = utils.plot() ax.set_ylim(0.0, 1.0) ax.plot(sample_times, sample_discounts); Market quotes 28 Also, here’s a bond priced by discounting its coupons on the curve: In : schedule = ql.Schedule(today, calendar.advance(today, 15, ql.Years), ql.Period(ql.Semiannual), calendar, ql.ModifiedFollowing, ql.ModifiedFollowing, ql.DateGeneration.Backward, False) bond = ql.FixedRateBond(3, 100.0, schedule, [0.04], ql.Actual360()) bond.setPricingEngine( ql.DiscountingBondEngine(ql.YieldTermStructureHandle(curve))) print(bond.cleanPrice()) Out: 105.77449622178025 “It looked like a good idea at the time” Now, let’s add an observer that checks whether the bond is out of date, and if so recalculates the bond and outputs its new price. In Python, I can do this by defining a function to be triggered by the notifications, by passing it to the observer I’m creating, and (this last step is as in C++) by registering the observer with the bond. As a reminder of how the whole thing works: the changes will come from the market quotes, but the observer doesn’t need to be concerned with that and only registers with the object it’s ultimately Market quotes 29 interested in; in this case, the bond whose price it wants to monitor. A change in any of the market quotes will cause the quote to notify the helper, which in turn will notify the curve, and so on to the pricing engine, the bond and finally our observer. In : prices = [] def print_price(): p = bond.cleanPrice() prices.append(p) print(p) o = ql.Observer(print_price) o.registerWith(bond) The function also appends the new price to a list that can be used later as a history of the prices. Let’s see if it works: In : quotes.setValue(101.0) Out: 105.77449622178025 105.86560413454633 Whoa, what was that? The function was called twice, which surprised me too when I wrote this notebook. It turns out that, due to a glitch of multiple inheritance, the curve sends two notifications to the instrument. After the first, the instrument recalculates but the curve doesn’t (which explains why the price doesn’t change); after the second, the curve updates and the price changes. This should be fixed in a future release, but again it doesn’t change the point of the notebook. Let’s set the quote back to its original value. In : quotes.setValue(100.0) Out: 105.86560413454633 105.77449632602422 Now, let’s say the market moves up and, accordingly, all the bonds prices increase to 101. Therefore, we need to update all the quotes. Market quotes 30 In : prices = [] for q in quotes: q.setValue(101.0) Out: 105.77449632602422 105.28388421801475 105.28388421801475 105.2186292238055 105.2186292238055 105.31959063332916 105.31959063332916 105.4878663261018 105.4878663261018 105.68032057288742 105.68032057288742 105.87580396665263 105.87580396665263 106.06201692071569 106.06201692071569 106.23044659597666 106.23044659597666 106.37409219514531 106.37409219514531 106.48708840535757 106.48708840535757 106.56505216839506 106.56505216839506 106.6057074004982 106.6057074004982 106.60980178159524 106.60980178159524 106.58011175261714 106.58011175261714 106.52070688437946 As you see, each of the updates sent a notification and thus triggered a recalculation. We can use the list of prices we collected (slicing it to skip duplicate values) to visualize how the price changed. In : unique_prices = prices[::2]+prices[-1::] _, ax = utils.plot() ax.plot(unique_prices, '-'); Market quotes 31 The first price is the original one, and the last price is the final one; but all those in between are calculated based on an incomplete set of changes in which some quotes were updated and some others weren’t. Those are all incorrect, and (since they went both above and below the range of the real prices) outright dangerous in case there were any triggers on price levels that could have fired. Clearly, this is not the kind of behavior we want our code to have. Alternatives? There are workarounds we can apply. For instance, it’s possible to freeze the bond temporarily, preventing it from forwarding notifications. In : bond.freeze() Now, notifications won’t be forwarded by the bond and thus won’t reach our observer. In fact, the following loop won’t print anything. In : for q in quotes: q.setValue(101.5) When we restore the bond, it sends a single notification, which triggers only one recalculation and gives the correct final price. Market quotes 32 In : bond.unfreeze() Out: 106.85839332472605 When using C++, it’s also possible to disable and re-enable notifications globally, which makes it more convenient. But it all feels a bit convoluted anyway. The whole thing will be simpler if we discard the initial idea and don’t force a recalculation for each notification. Pull, don’t push It’s preferable for updates to not trigger recalculation and just set some kind of dirty flag, just like the instruments in the library do. This way, you can control when the calculation occur. To do so, let’s remove the observer we have in place… In : del o …and instead create one that raises a flag when it’s notified. In : flag = {} flag['status'] = 'down' def set_flag(): flag['status'] = 'up' o = ql.Observer(set_flag) o.registerWith(bond) The flag is initially down… In : print(flag) Out: {'status': 'down'} …and quote changes cause it to be raised. In : for q in quotes: q.setValue(100.0) In : print(flag) Out: {'status': 'up'} At this point, we can ask the bond for its final price. Market quotes 33 In : bond.cleanPrice() Out: 105.77449633968729 Better yet, we can let the instrument do that: let’s remove the second observer, too, and just ask the instrument for its price after the changes. The instrument keeps track of whether it needs recalculation, so it doesn’t need us to keep track of it. In : del o In : for q in quotes: q.setValue(101.0) In : bond.cleanPrice() Out: 106.52070679688055 So, less is more? In this case, it seems so. 5. Term structures and their reference dates In this notebook, I show briefly how to set up term structures so that they track (or don’t track) the global evaluation date. Setup Import the QuantLib module and set up the global evaluation date. You might want to take note of the date, since we’ll be moving it around later on. In : import QuantLib as ql In : ql.Settings.instance().evaluationDate = ql.Date(3, ql.October, 2014) Specifying the reference date of a term structure In not-too-accurate terms, the reference date of a term structure is where it begins. It can be the evaluation date, but you might also want it to start on the spot date, for instance. We have two possibilities to define a reference date for a curve—even though some particular classes only allow one of them. The first is to define it by means of a (possibly null) offset from the current evaluation date; e.g., “two business days after the evaluation date” to define it as the spot date, or “no business days” to define it as the evaluation date itself. I’ll do it here by building a sample curve over a few swaps. Never mind the helper object that I’m building here… In : helpers = [ ql.SwapRateHelper(ql.QuoteHandle(ql.SimpleQuote(rate/100.0)), ql.Period(*tenor), ql.TARGET(), ql.Annual, ql.Unadjusted, ql.Thirty360(ql.Thirty360.BondBasis), ql.Euribor6M()) for tenor, rate in [((2,ql.Years), 0.201), ((3,ql.Years), 0.258), ((5,ql.Years), 0.464), ((10,ql.Years), 1.151), ((15,ql.Years), 1.588)] ] …because the construction of the curve is the main point: note the 0 and TARGET() arguments, specifying the number of days and the calendar used to determine business days. Term structures and their reference dates 35 In : curve1 = ql.PiecewiseFlatForward(0, ql.TARGET(), helpers, ql.Actual360()) The second possibility is to specify the reference date explicitly. For instance, the ForwardCurve class takes a vector of specific dates and the corresponding rates and interpolates between them; the first passed date is taken as the reference date of the curve. For comparison purposes, I’ll ask the curve above for its nodes and use them to build a ForwardCurve instance: In : dates, rates = zip(*curve1.nodes()) In : curve1.nodes() Out: ((Date(3,10,2014), 0.0019777694879293093), (Date(7,10,2016), 0.0019777694879293093), (Date(9,10,2017), 0.0036475517704509294), (Date(7,10,2019), 0.007660760701876805), (Date(7,10,2024), 0.018414773669420893), (Date(8,10,2029), 0.025311634328221498)) The curve built based on these data will be the same as the first, except that we’re specifying its reference date explicitly as October 3rd (the first passed date). In : curve2 = ql.ForwardCurve(dates, rates, ql.Actual360()) Both curves are defined over the same range of dates… In : print("{0} to {1}".format(curve1.referenceDate(), curve1.maxDate())) print("{0} to {1}".format(curve2.referenceDate(), curve2.maxDate())) Out: October 3rd, 2014 to October 8th, 2029 October 3rd, 2014 to October 8th, 2029 …and return the same rates, whether we ask for a given time (for instance, 5 years)… In : print(curve1.zeroRate(5.0, ql.Continuous)) print(curve2.zeroRate(5.0, ql.Continuous)) Out: 0.452196 % Actual/360 continuous compounding 0.452196 % Actual/360 continuous compounding …or for a given date. Term structures and their reference dates 36 In : print(curve1.zeroRate(ql.Date(7, ql.September, 2019), ql.Actual360(), ql.Continuous)) print(curve2.zeroRate(ql.Date(7, ql.September, 2019), ql.Actual360(), ql.Continuous)) Out: 0.452196 % Actual/360 continuous compounding 0.452196 % Actual/360 continuous compounding With the help of a couple more Python modules, we can also plot the whole curve by asking for rates over a set of times: In : %matplotlib inline import utils from matplotlib.ticker import FuncFormatter import numpy as np In : times = np.linspace(0.0, 15.0, 400) rates = [ curve1.zeroRate(t, ql.Continuous).rate() for t in times ] _, ax = utils.plot() ax.yaxis.set_major_formatter( FuncFormatter(lambda r,pos: utils.format_rate(r,2))) ax.plot(times, rates); Term structures and their reference dates 37 Moving the evaluation date To recap: we built the first curve specifying its reference date relative to the evaluation date, and the second curve specifying its reference date explicitly. Now, what happens if we change the evaluation date? In : ql.Settings.instance().evaluationDate = ql.Date(19, ql.September, 2014) As you might expect, the reference date of the first curve changes accordingly while that of the second curve doesn’t. We can see how the range of definition has now changed for the first curve, but not for the second: In : print("{0} to {1}".format(curve1.referenceDate(), curve1.maxDate())) print("{0} to {1}".format(curve2.referenceDate(), curve2.maxDate())) Out: September 19th, 2014 to September 24th, 2029 October 3rd, 2014 to October 8th, 2029 And of course the rates have changed, too… In : print(curve1.zeroRate(5.0, ql.Continuous)) print(curve2.zeroRate(5.0, ql.Continuous)) Out: 0.452196 % Actual/360 continuous compounding 0.452196 % Actual/360 continuous compounding …if we look at them in the right way. The whole curve has moved back a couple of weeks, so if we ask for a given time we’ll get the same rates; in other words, we’re asking for the zero rate over five years after the reference date, and that remains the same for a rigid translation of the curve. If we ask for the zero rate at a given date, though, we’ll see the effect: In : print(curve1.zeroRate(ql.Date(7, ql.September, 2019), ql.Actual360(), ql.Continuous)) print(curve2.zeroRate(ql.Date(7, ql.September, 2019), ql.Actual360(), ql.Continuous)) Out: 0.454618 % Actual/360 continuous compounding 0.452196 % Actual/360 continuous compounding Notifications Finally, we can see how the two curves behave differently also with respect to notifications. Let’s make two observers… Term structures and their reference dates 38 In : def make_observer(i): def say(): s = "Observer %d notified" % i print('-'*len(s)) print(s) print('-'*len(s)) return ql.Observer(say) obs1 = make_observer(1) obs2 = make_observer(2) …and check that they work correctly by connecting them to a few quotes. The first observer will receive notifications from the first and third quote, and the second observer will receive notifications from the second and third quote. In : q1 = ql.SimpleQuote(1.0) obs1.registerWith(q1) q2 = ql.SimpleQuote(2.0) obs2.registerWith(q2) q3 = ql.SimpleQuote(3.0) obs1.registerWith(q3) obs2.registerWith(q3) If I trigger a change in the first quote, the first observer is notified and outputs a message: In : q1.setValue(1.5) Out: ------------------- Observer 1 notified ------------------- A change in the second quote causes a message from the second observer… In : q2.setValue(1.9) Out: ------------------- Observer 2 notified ------------------- …and a change in the third quote causes both observers to react. Term structures and their reference dates 39 In : q3.setValue(3.1) Out: ------------------- Observer 1 notified ------------------- ------------------- Observer 2 notified ------------------- Now let’s connect the observers to the curves. The first observer will receive notifications from the curve that moves with the evaluation date, and the second observer will receive notifications from the curve that doesn’t move. In : obs1.registerWith(curve1) obs2.registerWith(curve2) Now we can see what happens when the evaluation date changes again: In : ql.Settings.instance().evaluationDate = ql.Date(23, ql.September, 2014) Out: ------------------- Observer 1 notified ------------------- As you can see, only the moving curve sent a notification. The other did not, since it was not modified by the change of evaluation date. 6. Pricing over a range of days Based on questions on Stack Exchange from Charles¹, bob.jonst², MCM³ and lcheng⁴. In : import QuantLib as ql import numpy as np np.random.seed(42) Let’s say we have an instrument (a fixed-rate bond, for instance) that we want to price on a number of dates. I assume we also have the market quotes, or the curves, corresponding to each of the dates; in this case we only need interest rates, but the library works the same way for any quotes. We’ll store the resulting prices in a dictionary, with the date as the key. In : prices = {} Producing a single price To price the bond on a single date, we create the instrument itself… In : start_date = ql.Date(8, ql.February, 2016) maturity_date = start_date + ql.Period(5, ql.Years) schedule = ql.Schedule(start_date, maturity_date, ql.Period(ql.Semiannual), ql.TARGET(), ql.Following, ql.Following, ql.DateGeneration.Backward, False) coupons = [0.01]*10 bond = ql.FixedRateBond(3, 100, schedule, coupons, ql.Thirty360(ql.Thirty360.BondBasis)) …and the required discount curve. For brevity, here I’m interpolating precomputed rates; I might as well bootstrap the curve on a set of market rates. ¹https://stackoverflow.com/questions/32869325/ ²https://quant.stackexchange.com/questions/35961/ ³https://quant.stackexchange.com/questions/38509 ⁴https://quant.stackexchange.com/questions/36830/ Pricing over a range of days 41 In : today = ql.Date(9, ql.May, 2018) nodes = [ today + ql.Period(i, ql.Years) for i in range(11) ] rates = [ 0.007, 0.010, 0.012, 0.013, 0.014, 0.016, 0.017, 0.018, 0.020, 0.021, 0.022 ] discount_curve = ql.ZeroCurve(nodes, rates, ql.Actual360()) Given the bond and the curve, we link them together through an engine, set the evaluation date and get the result. In : discount_handle = ql.RelinkableYieldTermStructureHandle(discount_curve) bond.setPricingEngine(ql.DiscountingBondEngine(discount_handle)) In : ql.Settings.instance().evaluationDate = today In : prices[today] = bond.cleanPrice() print(prices[today]) Out: 99.18942082987543 Pricing on multiple days We could repeat the above for all dates, but it goes against the grain of the library. The architecture (see chapter 2 of Implementing QuantLib⁵ for details) was designed so that the instrument can react to changing market conditions; therefore, we can avoid recreating the instrument. We’ll only change the discount curve and the evaluation date. For instance, here I’ll calculate the price for the business day before today: In : calendar = ql.TARGET() yesterday = calendar.advance(today, -1, ql.Days) I’ll generate random rates to avoid coming up with a new set; but the idea is to build the correct discount curve for the evaluation date. In : nodes = [ yesterday + ql.Period(i, ql.Years) for i in range(11) ] base_rates = np.array(rates) rates = base_rates * np.random.normal(loc=1.0, scale=0.005, size=base_rates.shape) discount_curve = ql.ZeroCurve(nodes, list(rates), ql.Actual360()) As I mentioned, I need to set the new evaluation date and to link the handle in the engine to the new discount curve… ⁵https://leanpub.com/implementingquantlib Pricing over a range of days 42 In : ql.Settings.instance().evaluationDate = yesterday discount_handle.linkTo(discount_curve) …after which the bond returns the updated price. In : prices[yesterday] = bond.cleanPrice() print(prices[yesterday]) Out: 99.16663635835845 By repeating the process, I can generate prices for, say, the whole of last year. Again, I’m generating random rates to avoid tedious listings or external data files; you’ll use the correct ones instead. In : first_date = calendar.advance(today, -1, ql.Years) date = calendar.advance(yesterday, -1, ql.Days) while date >= first_date: nodes = [ date + ql.Period(i, ql.Years) for i in range(11) ] rates = base_rates * np.random.normal(loc=1.0, scale=0.005, size=base_rates.shape) discount_curve = ql.ZeroCurve(nodes, list(rates), ql.Actual360()) ql.Settings.instance().evaluationDate = date discount_handle.linkTo(discount_curve) prices[date] = bond.cleanPrice() date = calendar.advance(date, -1, ql.Days) Here are the results. Through the random noise, you can see how the price increases as the bond gets nearer to maturity. In : %matplotlib inline import utils In : dates, values = zip(*sorted(prices.items())) In : fig, ax = utils.plot() ax.xaxis.set_major_formatter(utils.date_formatter()) ax.plot_date([ utils.to_datetime(d) for d in dates ], values,'-'); Pricing over a range of days 43 Using quotes If we work with quotes, we can also avoid rebuilding the curve. Let’s say our discount curve is defined as a risk-free curve with an additional credit spread. The risk-free curve is bootstrapped from a number of market rates; for simplicity, here I’ll use a set of overnight interest-rate swaps, but you’ll use whatever makes sense in your case. In : index = ql.Eonia() tenors = [ ql.Period(i, ql.Years) for i in range(1,11) ] rates = [ 0.010, 0.012, 0.013, 0.014, 0.016, 0.017, 0.018, 0.020, 0.021, 0.022 ] quotes = [] helpers = [] for tenor, rate in zip(tenors, rates): q = ql.SimpleQuote(rate) h = ql.OISRateHelper(2, tenor, ql.QuoteHandle(q), index) quotes.append(q) helpers.append(h) One thing to note: I’ll setup the curve so that it moves with the evaluation date. This means that I won’t pass an explicit reference date, but a number of business days and a calendar. Passing 0, as Pricing over a range of days 44 in this case, will cause the reference date of the curve to equal the evaluation date; passing 2, for instance, would cause it to equal the corresponding spot date. In : risk_free_curve = ql.PiecewiseFlatForward(0, ql.TARGET(), helpers, ql.Actual360()) Finally, I’ll manage credit as an additional spread over the curve: In : spread = ql.SimpleQuote(0.01) discount_curve = ql.ZeroSpreadedTermStructure( ql.YieldTermStructureHandle(risk_free_curve), ql.QuoteHandle(spread)) Now we can recalculate today’s price… In : prices = {} ql.Settings.instance().evaluationDate = today discount_handle.linkTo(discount_curve) prices[today] = bond.cleanPrice() print(prices[today]) Out: 96.50362161659807 …and as before, we go back; except this time we don’t need to build a new curve. Instead, we can set new values to the quotes and they will trigger the necessary recalculations. In : date = calendar.advance(today, -1, ql.Days) base_rates = np.array(rates) while date >= first_date: rates = base_rates * np.random.normal(loc=1.0, scale=0.005, size=base_rates.shape) for q, r in zip(quotes, rates): q.setValue(r) spread.setValue(spread.value()*np.random.normal(loc=1.0, scale=0.005)) ql.Settings.instance().evaluationDate = date prices[date] = bond.cleanPrice() date = calendar.advance(date, -1, ql.Days) Note that we didn’t create any new object in the loop; we’re only settings new values to the quotes. Again, here are the results: Pricing over a range of days 45 In : dates, values = zip(*sorted(prices.items())) fig, ax = utils.plot() ax.xaxis.set_major_formatter(utils.date_formatter()) ax.plot_date([ utils.to_datetime(d) for d in dates ], values,'-'); A complication: past fixings For instruments that depend on the floating rate, we might need some past fixings. This is not necessarily related to pricing on a range of dates: even on today’s date, we need the fixing for the current coupon. Let’s set the instrument up… Pricing over a range of days 46 In : forecast_handle = ql.YieldTermStructureHandle(risk_free_curve) index = ql.Euribor6M(forecast_handle) bond = ql.FloatingRateBond( 3, 100, schedule, index, ql.Thirty360(ql.Thirty360.BondBasis) ) bond.setPricingEngine(ql.DiscountingBondEngine(discount_handle)) In : ql.Settings.instance().evaluationDate = today for q, r in zip(quotes, base_rates): q.setValue(r) spread.setValue(0.01) …and try to price it. No joy. In : print(bond.cleanPrice()) Out: --------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) Cell In , line 1 ----> 1 print(bond.cleanPrice()) File /usr/local/lib/python3.10/site-packages/QuantLib/QuantLib.py:17938, in Bond.clea\ nPrice(self, *args) 17937 def cleanPrice(self, *args): > 17938 return _QuantLib.Bond_cleanPrice(self, *args) RuntimeError: Missing Euribor6M Actual/360 fixing for February 6th, 2018 Being in the past, the fixing can’t be retrieved from the curve. We have to store it into the index, after which the calculation works: In : index.addFixing(ql.Date(6, ql.February,2018), 0.005) print(bond.cleanPrice()) Out: 97.11939323923686 When pricing on a range of dates, though, we need to take into account the fact that the current coupon changes as we go back in time. These two dates will work… Pricing over a range of days 47 In : ql.Settings.instance().evaluationDate = ql.Date(1, ql.March, 2018) print(bond.cleanPrice()) ql.Settings.instance().evaluationDate = ql.Date(15, ql.February, 2018) print(bond.cleanPrice()) Out: 96.84331874622794 96.79054303973298 …but this one causes the previous coupon to be evaluated, and that requires a new fixing: In : ql.Settings.instance().evaluationDate = ql.Date(1, ql.February, 2018) print(bond.cleanPrice()) Out: --------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) Cell In , line 2 1 ql.Settings.instance().evaluationDate = ql.Date(1, ql.February, 2018) ----> 2 print(bond.cleanPrice()) File /usr/local/lib/python3.10/site-packages/QuantLib/QuantLib.py:17938, in Bond.clea\ nPrice(self, *args) 17937 def cleanPrice(self, *args): > 17938 return _QuantLib.Bond_cleanPrice(self, *args) RuntimeError: Missing Euribor6M Actual/360 fixing for August 4th, 2017 Once we add it, the calculation works again. In : index.addFixing(ql.Date(4, ql.August, 2017), 0.004) print(bond.cleanPrice()) Out: 96.98060241422583 (If you’re wondering how the calculation worked before, since this coupon belonged to the bond: on the other evaluation dates, this coupon was expired and the engine could skip it without needing to calculate its amount. Thus, its fixing didn’t need to be retrieved.) More complications: future prices What if we go forward in time, instead of pricing on past dates? For one thing, we’ll need to forecast curves in some way. One way is to imply them from today’s curves: I talk about implied curves in another notebook, so I won’t repeat myself here. Let’s assume we have implied rates and we can set them. Once we do, we can price in the future just as easily as we do in the past. As I write this, it’s May 19th 2018, and June 1st is in the future: Pricing over a range of days 48 In : ql.Settings.instance().evaluationDate = ql.Date(1, ql.June, 2018) print(bond.cleanPrice()) Out: 97.2126812565699 However, there’s another problem, as pointed out by Mariano Zeron⁶ in a post to the QuantLib mailing list. If we go further in the future, the bond will require—so to speak—future past fixings. In : ql.Settings.instance().evaluationDate = ql.Date(1, ql.June, 2019) print(bond.cleanPrice()) Out: --------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) Cell In , line 3 1 ql.Settings.instance().evaluationDate = ql.Date(1, ql.June, 2019) ----> 3 print(bond.cleanPrice()) File /usr/local/lib/python3.10/site-packages/QuantLib/QuantLib.py:17938, in Bond.clea\ nPrice(self, *args) 17937 def cleanPrice(self, *args): > 17938 return _QuantLib.Bond_cleanPrice(self, *args) RuntimeError: Missing Euribor6M Actual/360 fixing for February 6th, 2019 Here the curve starts on June 1st 2019, and cannot retrieve the fixing at the start of the corresponding coupon. One way out of this might be to forecast fixings off the current curve and store them: In : ql.Settings.instance().evaluationDate = ql.Date(1, ql.June, 2018) future_fixing = index.fixing(ql.Date(6,ql.February,2019)) print(future_fixing) index.addFixing(ql.Date(6,ql.February,2019), future_fixing) Out: 0.011387399107860378 This way, they will be retrieved in the same way as real past fixings. ⁶https://sourceforge.net/p/quantlib/mailman/message/35270917/ Pricing over a range of days 49 In : ql.Settings.instance().evaluationDate = ql.Date(1, ql.June, 2019) print(bond.cleanPrice()) Out: 98.30830224923507 Of course, you might forecast them in a better way: that’s up to you. And if you’re worried that this might interfere with pricing on today’s date, don’t: stored fixings are only used if they’re in the past with respect to the evaluation date. The fixing I’m storing below for February 3rd 2021 will be retrieved if the evaluation date is later… In : index.addFixing(ql.Date(3,ql.February,2021), 0.02) ql.Settings.instance().evaluationDate = ql.Date(1, ql.June, 2021) print(index.fixing(ql.Date(3,ql.February,2021))) Out: 0.02 …but it will be forecast from the curve when it’s after the evaluation date: In : ql.Settings.instance().evaluationDate = ql.Date(1, ql.June, 2020) print(index.fixing(ql.Date(3,ql.February,2021))) Out: 0.011367299732914539 7. A note on random numbers and dimensionality Setup Import QuantLib and the graphing module. In : %matplotlib inline import matplotlib.pyplot as plt import QuantLib as ql Also, define a helper function to make the notebook less verbose. In : def set_unit_square(ax): ax.axis('scaled') ax.set_xlim([0,1]) ax.set_ylim([0,1]) Covering a unit square Let’s say we want to extract points inside a unit square; that is, pairs of points in the domain (0, 1) × (0, 1). The dimensionality of the problem is 2, since we need 2 numbers, x and y , per each sample. With pseudo-random numbers, it doesn’t matter much: we can just extract single numbers and form pairs from them. In : rng = ql.MersenneTwisterUniformRng(42) In : xs = [] ys = [] for i in range(2047): xs.append(rng.next().value()) ys.append(rng.next().value()) In : fig = plt.figure(figsize=(8,8)) ax = fig.add_subplot(1,1,1) set_unit_square(ax) ax.plot(xs,ys,'o'); A note on random numbers and dimensionality 51 The same doesn’t hold for quasi-random numbers, for which each sample is correlated to the one that follows it in order to cover the domain evenly. We can see this by plotting the sequence of Sobol numbers generated to cover the 1-dimensional unit interval: A note on random numbers and dimensionality 52 In : fig = plt.figure(figsize=(12,4)) for i, n in enumerate([0,1,2,3, 4,5,6,7, 9,11,13,15, 19,23,27,31]): rng = ql.SobolRsg(1) xs = [ rng.nextSequence().value() for j in range(n) ] ax = fig.add_subplot(4, 4, i+1) ax.axis('scaled') ax.set_xlim([0,1]) ax.set_ylim([-0.1,0.1]) ax.set_xticks([]) ax.set_yticks([]) ax.plot(xs,*len(xs),'o') ax.text(0.0, 0.15, 'n = %d' % n) The points are not added randomly at all, but in a predetermined sequence. This ruins the random properties of the sequence when used with the wrong dimensionality. (You can also see how an even coverage is only obtained for a number of samples of the form n = 2i − 1 for some i.) In : rng = ql.SobolRsg(1) In : xs = [] ys = [] for i in range(2047): xs.append(rng.nextSequence().value()) ys.append(rng.nextSequence().value()) In : fig = plt.figure(figsize=(8,8)) ax = fig.add_subplot(1,1,1) set_unit_square(ax) ax.plot(xs,ys,'o'); A note on random numbers and dimensionality 53 To cover the domain correctly, we have to use the right dimensionality. A note on random numbers and dimensionality 54 In : rng = ql.SobolRsg(2) In : xs = [] ys = [] for i in range(2047): x,y = rng.nextSequence().value() xs.append(x) ys.append(y) In : fig = plt.figure(figsize=(8,8)) ax = fig.add_subplot(1,1,1) set_unit_square(ax) ax.plot(xs,ys,'o'); A note on random numbers and dimensionality 55 The pattern above covers the square evenly, and also causes projections on the two axes to have good coverage (which wouldn’t happen, for instance, with a regular placement in rows and columns; the projections of most points would coincide). It is also interesting to see how the coverage is built as the number of samples increase: A note on random numbers and dimensionality 56 In : fig = plt.figure(figsize=(12,9)) for i, n in enumerate([0,1,2,3, 4,5,6,7, 15,31,63,127]): rng = ql.SobolRsg(2) ax = fig.add_subplot(3, 4, i+1) ax.set_xticks([]) ax.set_yticks([]) if n == 0: continue points = [ rng.nextSequence().value() for j in range(n) ] xs,ys = zip(*points) ax.axis('scaled') ax.set_xlim([0,1]) ax.set_ylim([0,1]) ax.plot(xs,ys,'o') ax.text(0.0, 1.05, 'n = %d' % n) A note on random numbers and dimensionality 57 Dimensionality of Monte Carlo simulations The classes in the QuantLib Monte Carlo framework will check the dimensionality of the generators they’re given and will warn you if it’s not correct. It’s still up to you to find the correct one, while writing your engines (for details on that, you can check chapter 6 of Implementing QuantLib¹). For instance, let’s say that you want to simulate three correlated stocks; and for sake of simplicity, let’s say they follow the Black-Scholes process. You’ll build a process for each of them… In : today = ql.Date(27,ql.January,2018) ql.Settings.instance().evaluationDate = today risk_free = ql.YieldTermStructureHandle( ql.FlatForward(today, 0.01, ql.Actual360())) processes = [ ql.BlackScholesProcess( ql.QuoteHandle(ql.SimpleQuote(S)), risk_free, ql.BlackVolTermStructureHandle( ql.BlackConstantVol(today, ql.TARGET(), sigma, ql.Actual360()))) for S, sigma in [(100, 0.20), ( 80, 0.25), (110, 0.18)] ] …and a single multi-dimensional process that correlates them. In this case, the resulting process has three random drivers. In : rho = [[1.0, 0.6, 0.8], [0.6, 1.0, 0.4], [0.8, 0.4, 1.0]] process = ql.StochasticProcessArray(processes, rho) print(process.factors()) Out: 3 Now, let’s say that we want to simulate paths over four steps, starting from today and ending one year from now. Each sample of the Monte Carlo simulation will need three random number for each step, for a total of 12 random numbers. This is the dimensionality of the problem; and, as I mentioned, the framework will check it and complain if it doesn’t match. (Please bear with me as I build the several classes needed for random-numbers generation. If find yourself doing this, you might want to write a helper function, like I do here.) ¹https://leanpub.com/implementingquantlib A note on random numbers and dimensionality 58 In : def rng(dimensionality): return ql.GaussianRandomSequenceGenerator( ql.UniformRandomSequenceGenerator( dimensionality, ql.UniformRandomGenerator(42))) times = [0.25, 0.50, 0.75, 1.0] generator = ql.GaussianMultiPathGenerator(process, times, rng(10)) Out: --------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) Cell In , line 8 2 return ql.GaussianRandomSequenceGenerator( 3 ql.UniformRandomSequenceGenerator( 4 dimensionality, 5 ql.UniformRandomGenerator(42))) 7 times = [0.25, 0.50, 0.75, 1.0] ----> 8 generator = ql.GaussianMultiPathGenerator(process, times, rng(10)) File /usr/local/lib/python3.10/site-packages/QuantLib/QuantLib.py:24679, in GaussianM\ ultiPathGenerator.__init__(self, *args) 24678 def __init__(self, *args): > 24679 _QuantLib.GaussianMultiPathGenerator_swiginit(self, _QuantLib.new_Gaussia\ nMultiPathGenerator(*args)) RuntimeError: dimension (10) is not equal to (3 * 4) the number of factors times the \ number of time steps As you might expect, the thing works with the correct dimensionality: In : generator = ql.GaussianMultiPathGenerator(process, times, rng(12)) In : sample = generator.next().value() In : fig = plt.figure(figsize=(12,6)) ax = fig.add_subplot(1,1,1) ts = [0.0] + times y_min = 80 y_max = 110 for i in range(3): p, = ax.plot(ts, sample[i], label='Stock %d' % (i+1)) ax.plot(ts, sample[i], 'o', color=p.get_color()) y_min = min(y_min, min(sample[i])) y_max = max(y_max, max(sample[i])) ax.set_xlim(0.0-0.02, 1.0+0.02) ax.set_xticks(ts) ax.set_ylim(y_min-2, y_max+2) ax.legend(loc='best'); A note on random numbers and dimensionality 59 Interest-rate curves 8. EONIA curve bootstrapping In the next notebooks, I’ll reproduce the results of the paper by F. M. Ametrano and M. Bianchetti, Everything You Always Wanted to Know About Multiple Interest Rate Curve Bootstrapping but Were Afraid to Ask (April 2, 2013). The paper is available at SSRN: http://ssrn.com/abstract=2219548. In : %matplotlib inline import math import utils In : import QuantLib as ql In : today = ql.Date(11, ql.December, 2012) ql.Settings.instance().evaluationDate = today First try We start by instantiating helpers for all the rates used in the bootstrapping process, as reported in figure 25 of the paper. The first three instruments are three 1-day deposit that give us discounting between today and the day after spot. They are modeled by three instances of the DepositRateHelper class with a tenor of 1 day and a number of fixing days going from 0 (for the deposit starting today) to 2 (for the deposit starting on the spot date). In : helpers = [ ql.DepositRateHelper(ql.QuoteHandle(ql.SimpleQuote(rate/100)), ql.Period(1,ql.Days), fixingDays, ql.TARGET(), ql.Following, False, ql.Actual360()) for rate, fixingDays in [(0.04, 0), (0.04, 1), (0.04, 2)] ] Then, we have a series of OIS quotes for the first month. They are modeled by instances of the OISRateHelper class with varying tenors. They also require an instance of the Eonia class, which doesn’t need a forecast curve and can be shared between the helpers. EONIA curve bootstrapping 62 In : eonia = ql.Eonia() In : helpers += [ ql.OISRateHelper(2, ql.Period(*tenor), ql.QuoteHandle(ql.SimpleQuote(rate/100)), eonia) for rate, tenor in [(0.070, (1,ql.Weeks)), (0.069, (2,ql.Weeks)), (0.078, (3,ql.Weeks)), (0.074, (1,ql.Months))] ] Next, five OIS forwards on ECB dates. For these, we need to instantiate the DatedOISRateHelper class and specify start and end dates explicitly. In : helpers += [ ql.DatedOISRateHelper(start_date, end_date, ql.QuoteHandle(ql.SimpleQuote(rate/100)), eonia) for rate, start_date, end_date in [ ( 0.046, ql.Date(16,ql.January,2013), ql.Date(13,ql.February,2013)), ( 0.016, ql.Date(13,ql.February,2013), ql.Date(13,ql.March,2013)), (-0.007, ql.Date(13,ql.March,2013), ql.Date(10,ql.April,2013)), (-0.013, ql.Date(10,ql.April,2013), ql.Date(8,ql.May,2013)), (-0.014, ql.Date(8,ql.May,2013), ql.Date(12,ql.June,2013))] ] Finally, we add OIS quotes up to 30 years. In : helpers += [ ql.OISRateHelper(2, ql.Period(*tenor), ql.QuoteHandle(ql.SimpleQuote(rate/100)), eonia) for rate, tenor in [(0.002, (15,ql.Months)), (0.008, (18,ql.Months)), (0.021, (21,ql.Months)), (0.036, (2,ql.Years)), (0.127, (3,ql.Years)), (0.274, (4,ql.Years)), (0.456, (5,ql.Years)), (0.647, (6,ql.Years)), (0.827, (7,ql.Years)), (0.996, (8,ql.Years)), (1.147, (9,ql.Years)), (1.280, (10,ql.Years)), (1.404, (11,ql.Years)), (1.516, (12,ql.Years)), (1.764, (15,ql.Years)), (1.939, (20,ql.Years)), (2.003, (25,ql.Years)), (2.038, (30,ql.Years))] ] The curve is an instance of PiecewiseLogCubicDiscount (corresponding to the PiecewiseYield- Curve class in C++; I won’t repeat the argument for this choice made in section 4.5 of the paper). We let the reference date of the curve move with the global evaluation date, by specifying it as 0 days after the latter on the TARGET calendar. The day counter chosen is not of much consequence, as it is only used internally to convert dates into times. Also, we enable extrapolation beyond the maturity of the last helper; that is mostly for convenience as we retrieve rates to plot the curve near its far end. EONIA curve bootstrapping 63 In : eonia_curve_c = ql.PiecewiseLogCubicDiscount(0, ql.TARGET(), helpers, ql.Actual365Fixed()) eonia_curve_c.enableExtrapolation() To compare the curve with the one shown in figure 26 of the paper, we can retrieve daily overnight rates over its first two years and plot them: In : today = eonia_curve_c.referenceDate() end = today + ql.Period(2,ql.Years) dates = [ ql.Date(serial) for serial in range(today.serialNumber(), end.serialNumber()+1) ] rates_c = [ eonia_curve_c.forwardRate(d, ql.TARGET().advance(d,1,ql.Days), ql.Actual360(), ql.Simple).rate() for d in dates ] In : _, ax = utils.plot() utils.highlight_x_axis(ax) utils.plot_curve(ax, dates, [(rates_c,'-')], format_rates=True) However, we still have work to do. Out plot above shows a rather large bump at the end of 2012 which is not present in the paper. To remove it, we need to model properly the turn-of-year effect. EONIA curve bootstrapping 64 Turn-of-year jumps As explained in section 4.8 of the paper, the turn-of-year effect is a jump in interest rates due to an increased demand for liquidity at the end of the year. The jump is embedded in any quoted rates that straddles the end of the year and must be treated separately; the YieldTermStructure class allows this by taking any number of jumps, modeled as additional discount factors, and applying them at the specified dates. Our problem is to estimate the size of the jump. To simplify analysis, we turn to flat forward rates instead of log-cubic discounts; thus, we instantiate a PiecewiseFlatForward curve (corresponding to PiecewiseYieldCurve in C++). In : eonia_curve_ff = ql.PiecewiseFlatForward(0, ql.TARGET(), helpers, ql.Actual365Fixed()) eonia_curve_ff.enableExtrapolation() To show the jump more clearly, I’ll restrict the plot to the first 6 months: In : end = today + ql.Period(6,ql.Months) dates = [ ql.Date(serial) for serial in range(today.serialNumber(), end.serialNumber()+1) ] rates_ff = [ eonia_curve_ff.forwardRate(d, ql.TARGET().advance(d,1,ql.Days), ql.Actual360(), ql.Simple).rate() for d in dates ] In : _, ax = utils.plot() utils.highlight_x_axis(ax) utils.plot_curve(ax, dates, [(rates_ff,'-')], format_rates=True) EONIA curve bootstrapping 65 As we see, the forward ending at the beginning of January 2013 is out of line. In order to estimate the jump, we need to estimate a “clean” forward that doesn’t include it. A possible estimate (although not the only one) can be obtained by interpolating the forwards around the one we want to replace. To do so, we extract the values of the forwards rates and their corresponding dates. In : nodes = list(eonia_curve_ff.nodes()) If we look at the first few nodes, we can clearly see that the seventh is out of line. In : nodes[:9] Out: [(Date(11,12,2012), 0.00040555533025081675), (Date(12,12,2012), 0.00040555533025081675), (Date(13,12,2012), 0.00040555533047721286), (Date(14,12,2012), 0.00040555533047721286), (Date(20,12,2012), 0.0007604110692568178), (Date(27,12,2012), 0.0006894305026004767), (Date(3,1,2013), 0.0009732981324671213), (Date(14,1,2013), 0.0006728161005748453), (Date(13,2,2013), 0.000466380545910482)] EONIA curve bootstrapping 66 To create a curve that doesn’t include the jump, we replace the relevant forward rate with a simple average of the ones that precede and follow… In : nodes = (nodes, (nodes+nodes)/2.0) nodes[:9] Out: [(Date(11,12,2012), 0.00040555533025081675), (Date(12,12,2012), 0.00040555533025081675), (Date(13,12,2012), 0.00040555533047721286), (Date(14,12,2012), 0.00040555533047721286), (Date(20,12,2012), 0.0007604110692568178), (Date(27,12,2012), 0.0006894305026004767), (Date(3,1,2013), 0.000681123301587661), (Date(14,1,2013), 0.0006728161005748453), (Date(13,2,2013), 0.000466380545910482)] …and instantiate a ForwardCurve with the modified nodes. In : temp_dates, temp_rates = zip(*nodes) temp_curve = ql.ForwardCurve(temp_dates, temp_rates, eonia_curve_ff.dayCounter()) For illustration, we can extract daily overnight nodes from the doctored curve and plot them alongside the old ones: In : temp_rates = [ temp_curve.forwardRate(d, ql.TARGET().advance(d,1,ql.Days), ql.Actual360(), ql.Simple).rate() for d in dates ] In : _, ax = utils.plot() utils.highlight_x_axis(ax) utils.plot_curve(ax, dates, [(temp_rates,'-'), (rates_ff,'--')], format_rates=True) EONIA curve bootstrapping 67 Now we can estimate the size of the jump. As the paper hints, it’s more an art than a science. I’ve been able to reproduce closely the results of the paper by extracting from the two curves the forward rate over the two weeks around the end of the year: In : d1 = ql.Date(31,ql.December,2012) - ql.Period(1,ql.Weeks) d2 = ql.Date(31,ql.December,2012) + ql.Period(1,ql.Weeks) In : F = eonia_curve_ff.forwardRate(d1, d2, ql.Actual360(), ql.Simple).rate() F_1 = temp_curve.forwardRate(d1, d2, ql.Actual360(), ql.Simple).rate() print(utils.format_rate(F,digits=3)) print(utils.format_rate(F_1,digits=3)) Out: 0.082 % 0.067 % We want to attribute the whole jump to the last day of the year, so we rescale it according to (F − F1 ) · t12 = J · tJ where t12 is the time between the two dates and tJ is the time between the start and end date of the end-of-year overnight deposit. This gives us a jump quite close to the value of 10.2 basis points reported in the paper. EONIA curve bootstrapping 68 In : t12 = eonia_curve_ff.dayCounter().yearFraction(d1,d2) t_j = eonia_curve_ff.dayCounter().yearFraction(ql.Date(31,ql.December,2012), ql.Date(2,ql.January,2013)) J = (F-F_1)*t12/t_j print(utils.format_rate(J,digits=3)) Out: 0.101 % As I mentioned previously, the jump can be added to the curve as a corresponding discount factor 1/(1 + J · tJ ) on the last day of the year. The information can be passed to the curve constructor, giving us a new instance: In : B = 1.0/(1.0+J*t_j) jumps = [ql.QuoteHandle(ql.SimpleQuote(B))] jump_dates = [ql.Date(31,ql.December,2012)] eonia_curve_j = ql.PiecewiseFlatForward(0, ql.TARGET(), helpers, ql.Actual365Fixed(), jumps, jump_dates) Retrieving daily overnight rates from the new curve and plotting them, we can see the jump quite clearly: In : rates_j = [ eonia_curve_j.forwardRate(d, ql.TARGET().advance(d,1,ql.Days), ql.Actual360(), ql.Simple).rate() for d in dates ] In : _, ax = utils.plot() utils.highlight_x_axis(ax) utils.plot_curve(ax, dates, [(rates_ff,'-'), (rates_j,'o')], format_rates=True) EONIA curve bootstrapping 69 We can now go back to log-cubic discounts and add the jump. In : eonia_curve = ql.PiecewiseLogCubicDiscount(0, ql.TARGET(), helpers, ql.Actual365Fixed(), jumps, jump_dates) e