Monday, December 15, 2014

Odds and ends: indexes

Greetings, everybody.

Ok, I'll make it quick.

First of all: all the slides from the Düsseldorf user meeting are already available on the QuantLib site. Go get 'em. If you want some introduction, my previous post might help.

Then: this week I'm posting some content from the book appendix, namely, an overview of indexes.

Finally: the holiday season is coming.

It's likely that I won't be posting next week, so let me wish a merry Christmas to those who, like my family, celebrate it, a happy holiday to those who celebrate something else, and some holiday cheer to those who don't (but hey, you're missing out, people). I'll see you here after the festivities.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

Odds and ends: indexes

Like other classes such as Instrument and TermStructure, the Index class is a pretty wide umbrella: it covers concepts such as interest-rate indexes, inflation indexes, stock indexes—you get the drift.

Needless to say, the modeled entities are diverse enough that the Index class has very little interface to call its own. As shown in the listing below, all its methods have to do with index fixings. The isValidFixingDate method tells us whether a fixing was (or will be made) on a given date; the fixingCalendar method returns the calendar used to determine the valid dates; and the fixing method retrieves a fixing for a past date or forecasts one for a future date. The remaining methods deal specifically with past fixings: the name method, which returns an identifier that must be unique for each index, is used to index (pun not intended) into a map of stored fixings; the addFixing method stores a fixing (or many, in other overloads not shown here); and the clearFixing method clears all stored fixings for the given index.

Listing: interface of the Index class.
    class Index : public Observable {
        virtual ~Index() {}
        virtual std::string name() const = 0;
        virtual Calendar fixingCalendar() const = 0;
        virtual bool isValidFixingDate(const Date& fixingDate) const = 0;
        virtual Real fixing(const Date& fixingDate,
                            bool forecastTodaysFixing = false) const = 0;
        virtual void addFixing(const Date& fixingDate,
                               Real fixing,
                               bool forceOverwrite = false);
        void clearFixings();

Why the map, and where is it in the Index class? Well, we started from the requirement that past fixings should be shared rather than per-instance; if one stored, say, the 6-months Euribor fixing for a date, we wanted the fixing to be visible to all instances of the same index, and not just the particular one whose addFixing method we called. (Note that by "instances of the same index" I mean here instances of the same specific index, not of the same class, which might group different indexes; for instance, USDLibor(3*Months) and USDLibor(6*Months) are not instances of the same index; two different USDLibor(3*Months) are.) This was done by defining and using an IndexManager singleton behind the curtains. Smelly? Sure, as all singletons. An alternative might have been to define static class variables in each derived class to store the fixings; but that would have forced us to duplicate them in each derived class with no real advantage (it would be as much against concurrency as the singleton). In any case, this is one of the things we'll have to rethink in the next big QuantLib revision, in which we'll have to tackle concurrency. If these past years are any indication, we might expect it around 2020. No, just kidding. Maybe.

Since the returned index fixings might change (either because their forecast values depend on other varying objects, or because a newly available fixing is added and replaces a forecast) the Index class inherits from Observable so that instruments can register with its instances and be notified of such changes.

At this time, Index doesn't inherit from Observer, although its derived classes do (not surprisingly, since forecast fixings will almost always depend on some observable market quote). This was not an explicit design choice, but rather an artifact of the evolution of the code and might change in future releases. However, even if we were to inherit Index from Observer, we would still be forced to have some code duplication in derived classes, for a reason which is probably worth describing in more detail.

I already mentioned that fixings can change for two reasons. One is that the index depends on other observables to forecast its fixings; in this case, it simply registers with them (this is done in each derived class, as each class has different observables). The other reason is that a new fixing might be made available, and that's more tricky to handle. The fixing is stored by a call to addFixing on a particular index instance, so it seems like no external notification would be necessary, and that the index can just call the notifyObservers method to notify its observers; but that's not the case. As I said, the fixings is shared; if we store today's 3-months Euribor fixing, it will be available to all instances of such index, and thus we want all of them to be aware of the change. Moreover, instruments and curves might have registered with any of those Index instances, so all of them must send in turn a notification.

The solution is to have all instances of the same index communicate by means of a shared object; namely, we used the same IndexManager singleton that stores all index fixings. As I said, IndexManager maps unique index tags to sets of fixings; also, by making the sets instances of the ObservableValue class, it provides the means to register and receive notification when one or more fixings are added for a specific tag (this class will be described in a later post. You don't need the details here).

All pieces are now in place. Upon construction, any Index instance will ask IndexManager for the shared observable corresponding to the tag returned by its name method. When we call addFixings on, say, some particular 6-months Euribor index, the fixing will be stored into IndexManager; the observable will send a notification to all 6-months Euribor indexes alive at that time; and all will be well with the world.

However, C++ still throws a small wrench in our gears. Given the above, it would be tempting to call
in the Index constructor and be done with it. However, it wouldn't work; for the reason that in the constructor of the base class, the call to the virtual method name wouldn't be polymorphic. (If you're not familiar with the darker corners of C++: when the constructor of a base class is executed, any data members defined in derived classes are not yet built. Since any behavior specific to the derived class is likely to depend on such yet-not-existing data, C++ bails out and uses the base-class implementation of any virtual method called in the base-class constructor body.) From here stems the code duplication I mentioned a few paragraphs earlier; in order to work, the above method call must be added to the constructor of each derived index class. The Index class itself doesn't have a constructor (apart from the default one that the compiler provides).

As an example of a concrete class derived from Index, the listing below sketches the InterestRateIndex class.

Listing: sketch of the InterestRateIndex class.
    class InterestRateIndex : public Index, public Observer {
        InterestRateIndex(const std::string& familyName,
                          const Period& tenor,
                          Natural settlementDays,
                          const Currency& currency,
                          const Calendar& fixingCalendar,
                          const DayCounter& dayCounter);
        : familyName_(familyName), tenor_(tenor), ... {

        std::string name() const;
        Calendar fixingCalendar() const;
        bool isValidFixingDate(const Date& fixingDate) const {
            return fixingCalendar().isBusinessDay(fixingDate);
        Rate fixing(const Date& fixingDate,
                    bool forecastTodaysFixing = false) const;
        void update() { notifyObservers(); }

        std::string familyName() const;
        Period tenor() const;
        ... // other inspectors

        Date fixingDate(const Date& valueDate) const;
        virtual Date valueDate(const Date& fixingDate) const;
        virtual Date maturityDate(const Date& valueDate) const = 0;
        virtual Rate forecastFixing(const Date& fixingDate) const = 0;
        std::string familyName_;
        Period tenor_;
        Natural fixingDays_;
        Calendar fixingCalendar_;
        Currency currency_;
        DayCounter dayCounter_;
As you might expect, such class defines a good deal of specific behavior besides that inherited from Index. To begin with, it inherits from Observer, too, since Index doesn't. The InterestRateIndex constructor takes the data needed to specify the index: a family name, as in "Euribor", common to different indexes of the same family such as, say, 3-months and 6-months Euribor; a tenor that specifies a particular index in the family; and additional information such as the number of settlement days, the index currency, the fixing calendar, and the day-count convention used for accrual.

The passed data are, of course, copied into the corresponding data members; then the index registers with a couple of observables. The first is the global evaluation date; this is needed because, as I'll explain shortly, there's a bit of date-specific behavior in the class that is triggered when an instance is asked for today's fixing. The second observable is the one which is contained inside IndexManager and provides notifications when new fixings are stored. We can identify this observable here: the InterestRateIndex class has all the information needed to determine the index, so it can implement the name method and call it. However, this also means that classes deriving from InterestRateIndex must not override name; since the overridden method would not be called in the body of this constructor, they would register with the wrong notifier. Unfortunately, this can't be enforced in C++, which doesn't have a keyword like final in Java or sealed in C#; but the alternative would be to require that all classes derived from InterestRateIndex register with IndexManager, which is equally not enforceable, probably more error-prone, and certainly less convenient.

The other methods defined in InterestRateIndex have different purposes. A few implement the required Index and Observer interfaces; the simplest are update, which simply forwards any notification, fixingCalendar, which returns a copy of the stored calendar instance, and isValidFixingDate, which checks the date against the fixing calendar.
    std::string InterestRateIndex::name() const {
        std::ostringstream out;
        out << familyName_;
        if (tenor_ == 1*Days) {
            if (fixingDays_==0) out << "ON";
            else if (fixingDays_==1) out << "TN";
            else if (fixingDays_==2) out << "SN";
            else out << io::short_period(tenor_);
        } else {
            out << io::short_period(tenor_);
        out << " " <<;
        return out.str();

    Rate InterestRateIndex::fixing(
                           const Date& d,
                           bool forecastTodaysFixing) const {
        QL_REQUIRE(isValidFixingDate(d), ...);
        Date today = Settings::instance().evaluationDate();
        if (d < today) {
            Rate pastFixing =
            QL_REQUIRE(pastFixing != Null<Real>(), ...);
            return pastFixing;
        if (d == today && !forecastTodaysFixing) {
            Rate pastFixing = ...;
            if (pastFixing != Null<Real>())
                return pastFixing;
        return forecastFixing(d);

    Date InterestRateIndex::valueDate(const Date& d) const {
        QL_REQUIRE(isValidFixingDate(d) ...);
        return fixingCalendar().advance(d, fixingDays_, Days);
The name method is a bit more complicated. It stitches together the family name, a short representation of the tenor, and the day-count convention to get an index name such as "Euribor 6M Act/360" or "USD Libor 3M Act/360"; special tenors such as overnight, tomorrow-next and spot-next are detected so that the corresponding acronyms are used.

The fixing method contains the most logic. First, the required fixing date is checked and an exception is raised if no fixing was supposed to take place on it. Then, the fixing date is checked against today's date. If the fixing was in the past, it must be among those stored in the IndexManager singleton; if not, an exception is raised since there's no way we can forecast a past fixing. If today's fixing was requested, the index first tries looking for the fixing in the IndexManager and returns it if found; otherwise, the fixing is not yet available. In this case, as well as for a fixing date in the future, the index forecasts the value of the fixing; this is done by calling the forecastFixing method, which is declared as purely virtual in this class and implemented in derived ones. The logic in the fixing method is also the reason why, as I mentioned, the index registers with the evaluation date; the behavior of the index depends on the value of today's date, so it need to be notified when it changes.

Finally, the InterestRateIndex class defines other methods that are not inherited from Index. Most of them are inspectors that return stored data such as the family name or the tenor; a few others deal with date calculations. The valueDate method takes a fixing date and returns the starting date for the instrument that underlies the rate (for instance, the deposit underlying a LIBOR, which for most currencies starts two business days after the fixing date); the maturityDate method takes a value date and returns the maturity of the underlying instrument (e.g., the maturity of the deposit); and the fixingDate method is the inverse of valueDate, taking a value date and returning the corresponding fixing date. Some of these methods are virtual, so that their behavior can be overridden; for instance, while the default behavior for valueDate is to advance the given number of fixing days on the given calendar, LIBOR index mandates first to advance on the London calendar, then to adjust the resulting date on the calendar corresponding to the index currency. For some reason, fixingDate is not virtual; this is probably an oversight that should be fixed in a future release.

Aside: how much generalization?

Some of the methods of the InterestRateIndex class were evidently designed with LIBOR in mind, since that was the first index of that kind implemented in the library. On the one hand, this makes the class less generic than one would like: for instance, if we were to decide that the 5-10 years swap-rate spread were to be considered an interest-rate index in its own right, we would be hard-pressed to fit it to the interface of the base class and its single tenor method. But on the other hand, it is seldom wise to generalize an interface without having a couple of examples of classes that should implement it; and a spread between two indexes (being just that; a spread, not an index) is probably not one such class.

Liked this post? Share it:

Tuesday, December 9, 2014

Report from the QuantLib user meeting in Düsseldorf

Welcome back.

This last Thursday and Friday I was in Düsseldorf for the second QuantLib user meeting.

Man, was it good.

It was two days of great talks, with a good mixture of technical and financial content (and by the way, the slides are being collected and some are already available on the QuantLib page for the meeting). Thanks go to IKB, that sponsored the event and provided the venue, and especially to Michael von der Driesch and Peter Caspers, who organized the meeting and kept it running smoothly. I'm looking forward to next year. Really, we should have more of these events around Europe.

I wasn't giving a talk this year, so I just relaxed and enjoyed the show. It started with the keynote of Ferdinando Ametrano (@Ferdinando1970) on open-source finance: QuantLib, OpenGamma and Bitcoin, with the latter being the success story and having the lion's share in the presentation. (If you want to read Ferdinando's papers on Bitcoin, they're here and here). One thing stuck with me among those Ferdinando said about QuantLib: it isn't nearly as used in academia as we would like. If you teach in a university, let me know if there's anything I can do to help with this.

The second talk was by Roland Lichters, which told us of how speed is paramount in the CVA application he and his team at Quaternion develop. They fine-tuned QuantLib until they could price a swap on a given scenario in 30 microseconds (which I thought quite impressive) with a full analysis in about 30 seconds, but Bermudan swaptions were much slower; a full analysis would take some 50 minutes. Taking inspiration from the Longstaff-Schwartz method for American Monte Carlo, they now run an initial simulation to calculate swaption prices as a regression function. This initial time investment allows them to save time during the analysis, bringing its time back to the order of seconds.

In the last talk of the morning, Matthias Groncki presented the work he and Michael von der Driesch did at IKB on ABS. They wanted to make these instruments easier to simulate and price, and they knew that QuantLib could help since it had the basic pieces of the puzzle. Based on those, they wrote an implementation of the instrument in C++ and exposed it to Python. This allowed them to create an easier interface through IPython notebooks, providing graphs as well as numerical results as an aid to understanding.

After lunch, we restarted with Dirk Eddelbuettel (@eddelbuettel).  Let me add a personal note here: Dirk has been helping with the project and packaging the Debian version of QuantLib almost since day one, and we never had the chance to meet in person before. I was delighted we could do it this time.
Dirk gave us a description of his initial work on exporting QuantLib to R, and how this ultimately led to the development of Rcpp which is now used by more than 300 CRAN packages (or about 5% of the total). Other cool tools shown in his talk: Shiny and RMarkdown, which can use R packages (including RQuantLib) to provide interactive web-based reports. One thing I noticed is that most of the problems the RQuantLib package seems to have are shared with the C++ library, or rather caused by it: singletons, no time component in dates, no simple example for textbook calculations such as Black-Scholes. Eventually, we'll have to do something about it.

Next, Sebastian Schlenkrich on the hot topic of this meeting: the modeling of the basic spread between Libor and funding curves (that is, either the OIS curve or possibly a cross-currency curve for deals collateralized in another currency). I won't go into the details, which are in his (and André Miemiec's) paper at SSRN. I'll note, though, that his talk was the one with the most lively discussion at the end, with Sebastian, Daniel Duffy, Ferdinando and myself going over what could be the right design for including his results into QuantLib. No clear winner emerged yet, but the interest in the topic makes me guess that we'll work on it.

The final talk of the day was from Peter Caspers, who described some of his several contributions to QuantLib during the past year. The list is impressive, covering no-arbitrage SABR; ZABR and SVI interpolations; TSR CMS coupon pricers; CMS-spread coupon pricers; Credit Risk+ model; Gaussian 1D models; simulated annealing optimizer; Runge-Kutta ODE solver; and dynamically created Mersenne-Twister RNGs. With a gentle metaphorical nudge in my ribs, Peter noticed that most of them are in the experimental folder. I agree; we need some kind of triage process to look at the contents of that folder and move them in the core library. I'll explore some ideas after next release.

And with the end of this talk (also immortalized by Dirk on Twitter) we were off to dinner.

And there was evening and there was morning.

The second day started with a talk by Daniel Duffy, who described the latest iteration of his PDE and FDM framework. He made a number of interesting points. He told how, in a mathematical domain such as this, functional programming and generics are more suited to the task than object-oriented programming; and how C++11/14 made this style of programming easier with type inference (auto), lambda functions, and generic function objects (std::function). Using loosely coupled functions, instead of class interfaces, made it possible to separate clearly the different layers of the framework (PDEs, FDM, interface, and so on). Daniel's code might get into QuantLib at some point (the main problem being that we're still supporting older compilers that don't implement C++11) but I've also been talking Daniel into releasing his code as a project on its own. We'll see.

Next, Paolo Mazzocchi, who's been working with Ferdinando on modeling basis spreads. In discount-based legacy systems, synthetic deposits are required for maturities shorter than the Libor tenor. Their approach was to parameterize the spread between Libor and funding curve in that range and fit the parameters to the existing market quotes. The details are in the slides (already available): expect a paper soon.

The last talk of Friday morning was by Bernd Lewerenz, on the subject of a new pricing engine he contributed for Asian options based on a model by Jan Večeř (and that, serendipitously, I merged in the main repository just days before the meeting). The engine makes it possible to efficiently price arithmetic-average Asian options, widely used in commodity markets. Also, Bernd told us he's been driving the thing by means of the QuantLib Perl bindings, which makes this their first documented use.

Starting the afternoon talks: Eric Ehlers and his Reposit project. It is a rewrite of his gensrc generator for the QuantLib Excel bindings, it is implemented as a SWIG module, has a web site, and is a work in progress. Eric's talk had the merit of sparking two developments: it attracted some interest in rekindling the QuantLib Calc addin, and made me think that this is probably the time to check SWIG's support for shared_ptr and see if we can simplify our SWIG interfaces.

Next: a joint presentation by Klaus Spanderen and Johannes Göttker-Schnetmann on their quest for a calibrated local volatility. They are pretty well along it (in fact, so much that they changed he title of their talk to be a lot more optimistic than when they proposed it). Not only the details are in their presentation (and I wouldn't do them justice, so go read the slides), but their code is on GitHub, too.

The final talk was from André Miemiec, talking of an open field of research: the arbitrage-free and consistent modeling of swaptions, CMS and CMS spreads. He's been building on QuantLib to price CMS spread options. A simple extension of the existing code gave him a correct implementation on the first try; but there's still work to be done on speed before it can be used in a production setting. It was nice to see the library help getting a working prototype quickly, though.

And with this, we got to the end. All the comments I heard were enthusiastic (and mine, too).  As I said above: definitely something that should be replicated somewhere else.

Thanks again to the sponsor, the organizers, the speakers and all those who attended. I hope I'll see you again next year.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

Liked this post? Share it:

Monday, November 17, 2014

Odds and ends: date calculations

Welcome back.

My blog schedule continues to be erratic, I know. But November is traditionally the month when writers get all motivated, so I might be working on some new book content (no, my book, not a novel). In the meantime, this post contains a few sections I had already published in the PDF version of the drafts, somewhat reworked to include some new information.

In other news, the QuantLib User Meeting 2014 is drawing near. I won't be giving a talk like I did last year, but I'm looking forward to be there anyway. I'll report on the talks in a future post.

Finally: as you might know if you are subscribed to the developers mailing list, it turns out that QuantLib 1.4 doesn't work with Clang 3.5 and the newly-released Boost 1.57. I'll be putting out a 1.4.1 release to fix the problem; as you read this, it might already be available from our download page. If not, try again in a day or two (but only if you're using Clang: if not, you don't need to upgrade).

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

Odds and ends: date calculations

Date calculations are among the basic tools of quantitative finance. As can be expected, QuantLib provides a number of facilities for this task; I briefly describe some of them in the following subsections.

Dates and periods

An instance of the Date class represents a specific day such as November 15th, 2014—today's date as I write this post. This class provides a number of methods for retrieving basic information such as the weekday, the day of the month, or the year; static information such as the minimum and maximum date allowed (at this time, January 1st, 1901 and December 31st, 2199, respectively) or whether or not a given year is a leap year; or other information such as a date's Excel-compatible serial number or whether or not a given date is the last date of the month. The complete list of available methods and their interface is documented in the reference manual. No time information is included (although we've been talking about this).

Capitalizing on C++ features, the Date class also overloads a number of operators so that date algebra can be written in a natural way; for example, one can write expressions such as ++d, which advances the date d by one day; d + 2, which yields the date two days after the given date; d2 - d1, which yields the number of days between the two dates; d - 3*Weeks, which yields the date three weeks before the given date (and incidentally, features a member of the available TimeUnit enumeration, the other members being Days, Months, and Years); or d1 < d2, which yields true if the first date is earlier than the second one. The algebra implemented in the Date class works on calendar days; neither bank holidays nor business-day conventions are taken into account.

The Period class models lengths of time such as two days, three weeks, or five years by storing a TimeUnit and an integer. It provides a limited algebra and a partial ordering. For the non mathematically inclined, this means that two Period instances might or might not be compared to see which is the shorter; while it is clear that, say, 11 months are less than one year, it is not possible to determine whether 60 days are more or less than two months without knowing which two months. When the comparison cannot be decided, an exception is thrown.

And of course, even when the comparison seems obvious, we managed to sneak in a few surprises. For instance, the comparison
    Period(7,Days) == Period(1,Weeks)
returns true. It seems correct, right? Hold that thought.


Holidays and business days are the domain of the Calendar class. Several derived classes exist which define holidays for a number of markets; the base class defines simple methods for determining whether or not a date corresponds to a holiday or a business day, as well as more complex ones for performing tasks such as adjusting a holiday to the nearest business day (where "nearest" can be defined according to a number of business-day conventions, listed in the BusinessDayConvention enumeration) or advancing a date by a given period or number of business days.

It might be interesting to see how the behavior of a calendar changes depending on the market it describes. One way would have been to store in the Calendar instance the list of holidays for the corresponding market; however, for maintainability we wanted to code the actual calendar rules (such as "the fourth Thursday in November" or "December 25th of every year") rather than enumerating the resulting dates for a couple of centuries. Another obvious way would have been to use polymorphism and the Template Method pattern; derived calendars would override the isBusinessDay method, from which all others could be implemented. This is fine, but it has the shortcoming that calendars would need to be passed and stored in shared_ptrs. The class is conceptually simple, though, and is used frequently enough that we wanted users to instantiate it and pass it around more easily—that is, without the added verbosity of dynamic allocation.

The final solution was the one shown in the listing below.
    class Calendar {
        class Impl {
            virtual ~Impl() {}
            virtual bool isBusinessDay(const Date&) const = 0;
        boost::shared_ptr<Impl> impl_;
        bool isBusinessDay(const Date& d) const {
            return impl_->isBusinessDay(d);
        bool isHoliday(const Date& d) const {
            return !isBusinessDay(d);
        Date adjust(const Date& d,
                    BusinessDayConvention c = Following) const {
            // uses isBusinessDay() plus some logic
        Date advance(const Date& d,
                     const Period& period,
                     BusinessDayConvention c = Following,
                     bool endOfMonth = false) const {
            // uses isBusinessDay() and possibly adjust()
        // more methods
It is a variation of the pimpl idiom, also reminiscent of the Strategy or Bridge patterns; these days, the cool kids might call it type erasure, too. Long story short: Calendar declares a polymorphic inner class Impl to which the implementation of the business-day rules is delegated and stores a pointer to one of its instances. The non-virtual isBusinessDay method of the Calendar class forwards to the corresponding method in Calendar::Impl; following somewhat the Template Method pattern, the other Calendar methods are also non-virtual and implemented (directly or indirectly) in terms of isBusinessDay. (The same technique is used in a number of other classes, such as DayCounter in the next section or Parameter from this post.)

Derived calendar classes can provide specialized behavior by defining an inner class derived from Calendar::Impl; their constructor will create a shared pointer to an Impl instance and store it in the impl_ data member of the base class. The resulting calendar can be safely copied by any class that need to store a Calendar instance; even when sliced, it will maintain the correct behavior thanks to the contained pointer to the polymorphic Impl class. Finally, we can note that instances of the same derived calendar class can share the same Impl instance. This can be seen as an implementation of the Flyweight pattern—bringing the grand total to about two and a half patterns for one deceptively simple class.

Enough with the implementation of Calendar, and back to its behavior. Here's the surprise I mentioned in the previous section. Remember Period(1,Weeks) being equal to Period(7,Days)? Except that for the advance method of a calendar, 7 days means 7 business days. Thus, we have a situation in which two periods p1 and p2 are equal (that is, p1 == p2 returns true) but calendar.advance(p1) differs from calendar.advance(p2).

Yay, us.

I'm not sure I have a good idea for a solution here. If we want backwards compatibility, the current uses of Days must keep working in the same way; so it's not possible, say, to start interpreting calendar.advance(7, Days) as 7 calendar days. One way out might be to keep the current situation, introduce two new enumeration cases BusinessDays and CalendarDays that remove the ambiguity, and deprecate Days. Another is to just remove the inconsistency by dictating that a 7-days period do not, in fact, equal one week; I'm not overly happy about this one.

If we give up on backwards compatibility (the legendary QuantLib 2.0) then there are more possibilities. One is to always use Days as calendar days and add BusinessDays as a different enumeration case. Another, which I'm liking more and more as I think about it, would be to always use Days as calendar days and add a specific method advanceBusinessDays to Calendar (or maybe an overload advance(n, BusinessDays), with BusinessDays being an instance of a separate class); however, this would mean that 3 business days wouldn't be a period.

As I said, no obvious solution. If you have any other suggestions, I'm all ears.

Day-count conventions

The DayCounter class provides the means to calculate the distance between two dates, either as a number of days or a fraction of an year, according to different conventions. Derived classes such as Actual360 or Thirty360 exist; they implement polymorphic behavior by means of the same technique used by the Calendar class and described in the previous section.

Unfortunately, the interface has a bit of a rough edge. Instead of just taking two dates, the yearFraction method is declared as
    Time yearFraction(const Date&,
                      const Date&,
                      const Date& refPeriodStart = Date(),
                      const Date& refPeriodEnd = Date()) const;
The two optional dates are required by one specific day-count convention (namely, the ISMA actual/actual convention) that requires a reference period to be specified besides the two input dates. To keep a common interface, we had to add the two additional dates to the signature of the method for all day counters (most of which happily ignore them). This is not the only mischief caused by this day counter; you'll see another in the next section.


The Schedule class, shown in the listing below, is used to generate sequences of coupon dates.
    class Schedule {
        Schedule(const Date& effectiveDate,
                 const Date& termination Date,
                 const Period& tenor,
                 const Calendar& calendar,
                 BusinessDayConvention convention,
                 BusinessDayConvention terminationDateConvention,
                 DateGeneration::Rule rule,
                 bool endOfMonth,
                 const Date& firstDate = Date(),
                 const Date& nextToLastDate = Date());
        Schedule(const std::vector<Date>&,
                 const Calendar& calendar = NullCalendar(),
                 BusinessDayConvention convention = Unadjusted);

        Size size() const;
        bool empty() const;
        const Date& operator[](Size i) const;
        const Date& at(Size i) const;
        const_iterator begin() const;
        const_iterator end() const;

        const Calendar& calendar() const;
        const Period& tenor() const;
        bool isRegular(Size i) const;
        Date previousDate(const Date& refDate) const;
        Date nextDate(const Date& refDate) const;
        ... // other inspectors and utilities
Following practice and ISDA conventions, it has to accept a lot of parameters; you can see them as the argument list of its constructor. (Oh, and you'll forgive me if I don't go and explain all of them. I'm sure you can guess what they mean.) They're probably too many, which is why the library uses the Named Parameter Idiom (already described in this post) to provide a less unwieldy factory class. With its help, a schedule can be instantiated as
    Schedule s = MakeSchedule().from(startDate).to(endDate)
Other methods include on the one hand, inspectors for the stored data; and on the other hand, methods to give the class a sequence interface, e.g., size, operator[], begin, and end.

The Schedule class has a second constructor, taking a precomputed vector of dates. It's only kind of working, though: the resulting Schedule instances simply don't have some of the data that their inspectors are supposed to return, so those methods will throw an exception if called. Among those, there are tenor and isRegular, about which I need to spend a couple of words.

First of all, isRegular(i) doesn't refer to the i-th date, but to the i-th interval; that is, the one between the i-th and (i+1)-th dates. This said, what does "regular" means? When a schedule is built based on a tenor, most intervals correspond to the passed tenor (and thus are regular) but the first and last intervals might be shorter or longer depending on whether we passed an explicit first or next-to-last date. We might do this, e.g., when we want to specify a short first coupon.

If we build the schedule with a precomputed set of dates, we don't have the tenor information and we can't tell if a given interval is regular. (Well, we could use heuristics, but it could get ugly fast.) The bad news is, this makes it impossible to use that schedule to build a sequence of bond coupons; if we pass it to the constructor of, say, a fixed-rate bond, we'll get an exception. And why, oh, why does the bond needs the missing info in order to build the coupons? Because the day-count convention of the bond might be ISMA actual/actual, which needs a reference period; and in order to calculate the reference period, we need to know the coupon tenor.

Fortunately, it shouldn't be difficult to fix this problem. One way would be to check the day-count convention, and only try to calculate the reference period when needed; this way the constructor would still raise an exception for ISMA actual/actual, but would succeed for all other conventions. Another way might be to add the tenor and regularity information to the schedule, so that the corresponding methods can work; but I'm not sure that this makes a lot of sense.

Liked this post? Share it:

Monday, October 13, 2014

Announcement: QuantLib user meeting 2014


There's two good news in this post.

The first is that, after preparing and recording about 6 hours of QuantLib content for Quants Hub, I'm back to regular sleep hours.

The video will be available for purchase from the Quants Hub site once it's been edited. I'll keep you posted.

The second is that, thanks to the sponsorship of IKB, we'll have a second QuantLib user group meeting in Dusseldorf on December 4th and 5th. The first was great: you can go and read about it in a couple of my old posts (here is a summary, and here is a screencast of my talk). Even better news: the registration is free and still open. More information is at this page on the QuantLib site. I'll be there, listening to the many excellent speakers.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

Liked this post? Share it:

Wednesday, October 1, 2014

QuantLib notebook: implied term structures

Hello everybody.

As I said in my previous post, I'm preparing for a workshop I'll be giving for Quants Hub (it's in London, October 8th, thanks for asking. Yes, you can still register.) This doesn't leave me any time for posting new book content, but since I'm rehearsing the material I thought I might as well record some of it. Thus, here is the third QuantLib notebook; this one is about implied term structures, and acknowledgments go to user Lisa Ann, who prompted it with her question.

In other news: remember last year's QuantLib User Meeting in Düsseldorf? We're having another in December. More details to come.

Ok, here's the video. Make it full screen, sit back, and enjoy.

(The video is also available on Vimeo.)

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

Liked this post? Share it:

Thursday, September 11, 2014

Still on hiatus, and some news

Welcome back, everybody. I hope you had a nice vacation.

I'm sorry, but it's going to be a while before I add new content to the blog. I'm preparing for a workshop I'll be giving for Quants Hub, and that's taking all the time I would usually give to the book and the rest. The workshop is going to be in London on October 8th. I hope to resume posting after that.

(For those interested: you can register for the workshop here. It's going to be different from the course I usually teach in collaboration with MoneyScience. On the one hand, it will not only focus on the development and extension of the library, but also on its usage. On the other hand, it will also be taped; the videos will be made available for purchase on the Quants Hub site at some date, which might come handy if, for instance, there's an ocean between you and London. I'll keep you informed.)

Another bit of update: the two videos that I've posted so far in the QuantLib Notebooks series are also available on YouTube now, so you can watch them there if you prefer it to Vimeo. You can also have a look at the series page on LeanPub if you didn't already.

Ok, enough advertising already. I hope to be back with some content soon.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

Liked this post? Share it:

Monday, July 14, 2014

Intermission: book update and summer schedule

Hello everybody, and congratulations to any Germans that might be reading. Nice job there in Rio.

No new content this week. Instead, I've updated the PDF version of my book with the changes I did while I was revising chapter 3, chapter 4 and part of appendix A to post them on the blog. Also, I'm thinking of making an ebook version available through LeanPub: if you're interested, visit the page I set up on their site and leave your feedback. There's a page for the QuantLib Notebooks, too.

Going onwards: my summer schedule for the blog will be, well, no fixed schedule. I might post if I find something interesting, but on the whole I think I'll take a step back and relax a bit. I'll try to be back in autumn with some new material.

On the other hand, the good Jacob Bettany and his crew at MoneyScience will be available the whole summer to accept your registrations for my Introduction to QuantLib Development course. It's three days in London (September 22nd, 23rd and 24th) in which I present material from Implementing QuantLib, you get to ask me questions and generally pick my brain, and I retaliate by making you do coding exercises and implement new instruments. I look forward to see you there.

I guess that's all for now, folks. Enjoy your summer holidays.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

Liked this post? Share it:

Monday, July 7, 2014

Chapter 6, part 8 of 8: example


This is the last post in a series (started here) that covered chapter 6 of my book, viz., the Monte Carlo framework in QuantLib. I've no idea yet what to write in next post: we'll see.

And of course, the weekly plug: you can still register for my Introduction to QuantLib Development course. It's three days of presentations (for me) and exercises (for you) based on the material in my book. It's going to be in London from September 22nd to 24th. Click the link above for more info.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

Example: basket option

To close this series, I'll show and discuss an example of how to build a pricing engine with the Monte Carlo machinery I described so far. The instrument I'll use is a simple European option on a basket of stocks, giving its owner the right to buy or sell the whole basket at an given price; the quantities of each of the stocks in the basket are also specified by the contract.

For brevity, I won't show the implementation of the instrument class itself, BasketOption. (This class is not the same as the BasketOption class implemented in QuantLib. The one used here is simplified for illustration purposes.) It would be quite similar to the VanillaOption class I've shown in this post, with the addition of a data member for the quantities (added to both the instrument and its arguments class). Also, I won't deal with the Greeks; but if the BasketOption class were to define them, methods such as delta and gamma would return a vector.

The main class in this example is the one implementing the pricing engine. It is the MCEuropeanBasketEngine class template, shown in listing 6.16.

Listing 6.16: Implementation of the McEuropeanBasketEngine class template.
    template <class RNG = PseudoRandom, class S = Statistics>
    class MCEuropeanBasketEngine
        : public BasketOption::engine,
          private McSimulation<MultiVariate,RNG,S> {
        typedef McSimulation<MultiVariate,RNG,S> simulation_type;
        typedef typename simulation_type::path_generator_type
        ... // same for the other defined types
                const shared_ptr<StochasticProcess>&,
                const Handle<YieldTermStructure>& discountCurve,
                Size timeSteps,
                Size timeStepsPerYear,
                bool antitheticVariate,
                Size requiredSamples,
                Real requiredTolerance,
                Size maxSamples,
                BigNatural seed);
        void calculate() const {
            const S& stats = this->mcModel_->sampleAccumulator();
            results_.value = stats.mean();
            if (RNG::allowsErrorEstimate)
                results_.errorEstimate = stats.errorEstimate();
        TimeGrid timeGrid() const;
        shared_ptr<path_generator_type> pathGenerator() const;
        shared_ptr<path_pricer_type> pathPricer() const;
        shared_ptr<StochasticProcess> process_;
        Handle<YieldTermStructure> discountCurve_;
        Size timeSteps_, timeStepsPerYear_;
        Size requiredSamples_;
        Size maxSamples_;
        Real requiredTolerance_;
        BigNatural seed_;

As expected, it inherits publicly from the BasketOption::engine class; however, it also inherits privately from the McSimulation class template.

In idiomatic C++, the use of private inheritance denotes an "is implemented in terms of" relationship. We don't want public inheritance here, since that would imply an "is a" relationship; in our conceptual model, MCEuropeanBasketEngine is a pricing engine for the basket option and not a simulation that could be used on its own. The use of multiple inheritance is shunned by some, and in fact it could be avoided in this case; our engine might use composition instead, and contain an instance of McSimulation. However, that would require inheriting a new simulation class from McSimulation in order to implement its purely virtual methods (such as pathGenerator or pathPricer) and would have the effect to make the design of the engine more complex; whereas, by inheriting from McSimulation, we can define such methods in the engine itself. (In order to avoid this dilemma, we would have to rewrite the McSimulation class so that it doesn't use the Template Method pattern; that is, we should make it a concrete class which would be passed the used path generator and pricer as constructor arguments. Apart from breaking backward compatibility, I'm not sure this would be worth the hassle.) Finally, note the template parameters: we leave to the user the choice of the RNG traits, but since we're modeling a basket option we specify the MultiVariate class as the Monte Carlo traits.

The constructor (whose implementation is not shown for brevity) takes the stochastic process for the underlyings, a discount curve, and a number of parameters related to the simulation. It copies the arguments into the corresponding data members (except for the antitheticVariate flag, which is passed to the McSimulation constructor) and registers the newly-built instance as an observer of both the stochastic process and the discount curve.

The calculate method, required by the PricingEngine interface, calls the method by the same name from McSimulation and passes it the needed parameters; then, it retrieves the statistics and stores the mean value and, if possible, the error estimate. This behavior is not specific of this particular option, and in fact it could be abstracted out in some generic McEngine class; but this would muddle the conceptual model. Such a class would inherit from McSimulation, but it could be reasonably expected to inherit from PricingEngine, too. However, it we did that, our engine would inherit from both McEngine and BasketOption::engine, leading to the dreaded inheritance diamond. On the other hand, if we didn't inherit it from PricingEngine (calling it McEngineAdapter or something) it would add more complexity to the inheritance hierarchy, since it would be an additional layer between McSimulation and our engine, and wouldn't remove all the scaffolding anyway: our engine would still need to define a calculate method forwarding to the one in the adapter. Thus, I guess we'll leave it at that.

The rest of listing 6.16 (continued below) shows the three methods required to implement the McSimulation interface and turn our class into a working engine.

    template <class RNG, class S>
    TimeGrid MCEuropeanBasketEngine<RNG,S>::timeGrid() const {
        Time T = process_->time(arguments_.exercise->lastDate());
        if (timeSteps_ != Null<Size>()) {
            return TimeGrid(T, timeSteps_);
        } else if (timeStepsPerYear_ != Null<Size>()) {
            Size steps = timeStepsPerYear_*T;
            return TimeGrid(T, std::max<Size>(steps, 1));
        } else {
            QL_FAIL("time steps not specified");
The timeGrid method builds the grid used by the simulation by specifying the end time (given by the exercise date) and the number of steps. Depending on the parameters passed to the engine constructor, they might have been specified as a given total number (the first if clause) or a number per year (the second). In either case, the specification is turned into a total number of steps and passed to the TimeGrid constructor. If neither was specified (the else clause) an error is raised.

    template <class RNG, class S>
    shared_ptr<typename MCEuropeanBasketEngine<RNG,S>::
    MCEuropeanBasketEngine<RNG,S>::pathGenerator() const {
        Size factors = process_->factors();
        TimeGrid grid = timeGrid();
        Size steps = grid.size() - 1;
        typename RNG::rsg_type gen =
            RNG::make_sequence_generator(factors*steps, seed_);
        return shared_ptr<path_generator_type>(
                   new path_generator_type(process_, grid, gen));
The pathGenerator method asks the process for the number of factors, determines the number of time steps from the result of the timeGrid method I just described, builds a random-sequence generator with the correct dimensionality (which of course is the product of the two numbers above) and uses it together with the underlying process and the time grid to instantiate a multi-path generator.

    template <class RNG, class S>
    shared_ptr<typename MCEuropeanBasketEngine<RNG,S>::
    MCEuropeanBasketEngine<RNG,S>::pathPricer() const {
        Date maturity = arguments_.exercise->lastDate();
        return shared_ptr<path_pricer_type>(
            new EuropeanBasketPathPricer(
                arguments_.payoff, arguments_.quantities,
Finally, the pathPricer method collects the data required to determine the payoff on each path—that is, the payoff object itself, the quantities and the discount at the exercise date, that the method precalculates—and builds an instance of the EuropeanBasketPathPricer class, sketched in listing 6.17.

Listing 6.17: Sketch of the EuropeanMultiPathPricer class.
    class EuropeanBasketPathPricer : public PathPricer<MultiPath> {
        EuropeanBasketPathPricer(const shared_ptr<Payoff>& payoff,
                                 const vector<Real>& quantities,
                                 DiscountFactor discount);
        Real operator()(const MultiPath& path) const {
            Real basketValue = 0.0;
            for (Size i=0; i<quantities_.size(); ++i)
                basketValue += quantities_[i] * path[i].back();
            return (*payoff)(basketValue) * discount_;
        shared_ptr<Payoff> payoff_;
        vector<Real> quantities_;
        DiscountFactor discount_;

The path pricer stores the passed data and uses them to calculate the realized value of the option in its operator() goverloading, whose implementation calculates the value of the basket at maturity and returns the corresponding payoff discounted to the present time. The calculation of the basket value is a straightforward loop that combines the stored quantities with the asset values at maturity, retrieved from the end points of the respective paths.

At this point, the engine is completed; but it's still a bit unwieldy to instantiate. We'd want to add a factory class with a fluent interface, so that one can write
    engine = MakeMcEuropeanBasketEngine<PseudoRandom>(process,
but I'm not showing its implementation here. There's plenty of such examples in the library for your perusal.

Also, you'll have noted that I kept the example as simple as possible, and for that reason I avoided a few possible generalizations. For instance, another pricer might need some other specific dates besides the maturity, and the time grid should be built accordingly; or the payoff might have been path-dependent, so that the path pricer should look at several path values besides the last. But I don't think you'll have any difficulties in writing the corresponding code.

Aside: need-to-know basis

As usual, I didn't show in which files one should put the several classes I described. Well, as a general principle, the more encapsulation the better; thus, to begin with, helper classes should be hidden from the public interface. For instance, the EuropeanBasketPathPricer class is only used by the engine and could be hidden from client code. The best way would be to define such classes inside an anonymous namespace in a .cpp file, but that's not always possible: in our case, the engine is a class template, so that's not an option. The convention we're using in QuantLib is that helper classes that must be declared in a header file are placed in a nested namespace detail, and the user agrees to be a gentleman (or a lady, of course) and to leave it alone.

If, later on, we find out that we need the class elsewhere (for instance, because we want to use EuropeanBasketPathPricer as a control-variate path pricer for another engine) we can move it to a header file—or, if it's already in one, to the main QuantLib namespace—and make it accessible to everybody.

After all this, you might still have a question: how generic is this engine, really? Well, somewhat less that I'd like. It is generic in the sense that you can plug in any process for N asset prices, and it will work. It can even do more exotic stuff, such as quanto effects: if you go the traditional way and model the effect as an correction for the drift of the assets, you can write (or better yet, decorate) your process so that its drift or evolve method takes into account the correction, and you'll be able to use it without changing the engine code.

However, if you want to use a process that models other stochastic variables besides the asset prices (say, a multi-asset Heston process) you're likely to get into some trouble. The problem is that you'll end up in the operator() of the path pricer with N quantities and 2N paths, without any indication of which are for the prices and which for the volatilities. How should they be combined to get the correct basket price? Of course, one can write a specific pricer for any particular process; but I'd like to provide something more reusable.

One simple solution is for the process and the pricer to share some coding convention. For instance, if you arrange the process so that the first N paths are those of the prices and the last M are those of any other variables, the pricer in listing 6.17 will work; note that it loops over the number of quantities, not the size of the process. However, this leaves one wide open to errors that will go undetected if a process doesn't conform to the convention.

Another possibility that comes to mind is to decorate the process with a layer that would show the asset prices and hide the extra variables. The decorator would appear as a process requiring the same number of random variates as the original, but exposing less stochastic variables. It would have to store some state in order to keep track of the hidden variables (which wouldn't be passed to its methods): its evolve method, to name one, would have to take the current prices, add the extra variables stored in the instance, call evolve in the original process, store the new values of the extra variables so that they're available for next call, and return only the new asset prices.

However, I'm not a fan of this solution. It would only work when the decorated methods are called in the expected order, and would break otherwise; for instance, if you called a method twice with the exact same arguments (starting prices and random variates) you would get different return, due to the changed internal state. The long way to express this is that the methods of the decorated process would lose referential transparency. The short way is that its interface would be a lie.

A more promising possibility (albeit one that requires some changes to the framework) would be to let the path generator do the filtering, with some help from the process. If the process could provide a method returning some kind of mask—say, a list of the indices of the assets into the whole set of variables, or a vector of booleans where the i-th element would be set to true if the i-th path corresponds to an asset price—then the generator could use it to drive the evolution of the variables correctly, and at the same time write in the multi-path only the asset prices. To make the whole thing backward-compatible, the mask would be passed to the generator only optionally, and the base StochasticProcess class would have a default implementation of the method returning an empty mask. If either way no mask were provided, the generator would revert to the current behavior. If needed, the method could be generalized to return other masks besides that for the prices; for instance, a variance swap would be interested in the volatility but not the price.

Finally, note that some assumptions are built right into the engine and cannot be relaxed without rewriting it. For instance, we're assuming that the discount factor for the exercise date is deterministic. If you wrote a process that models stochastic interest rates as well as asset prices, and therefore wanted to discount on each path according to the realization of the rates on that path, you'd have to write a custom engine and path pricer to be used with your specific process. Fortunately, that would be the only things you'd have to write; you'd still be able to reuse the other pieces of the framework.

Liked this post? Share it:

Monday, June 30, 2014

Chapter 6, part 7 of 8: Monte Carlo simulations

Hello everybody.

This is the seventh in a series of posts that started here and covers chapter 6 of my book; that is, the Monte Carlo framework.

A week or two ago, I've been issuing a call to arms on the QuantLib mailing list. We have quite a few old items in our bug tracker and patch tracker, and we should see if they're still relevant. A few have been closed already, but if you want to get in the action and help triaging the rest, please write to the mailing list.

Do I still need to tell you about my Introduction to QuantLib Development course? Apparently I do. It's three days of lectures and exercises based on Implementing QuantLib, it's in London from September 22nd to 24th, and more information, a brochure and registration form are available here.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the widgets for that are in the sidebar, at the top right of the page. Also, make sure to check my Training page.

Monte Carlo simulations

Up one step in complexity from the Monte Carlo model of last post and we get to the McSimulation class template, sketched in listing 6.15. Its job is to drive the simple-minded MonteCarloModel towards a goal, be it a required accuracy or number of samples. It can be used (and it was designed) as a starting point to build a pricing engine.

Listing 6.15: Sketch of the McSimulation class template.
    template <template <class> class MC, class RNG,
              class S = Statistics>
    class McSimulation {
        typedef typename MonteCarloModel<MC,RNG,S>::path_generator_type
        typedef typename MonteCarloModel<MC,RNG,S>::path_pricer_type
        typedef typename MonteCarloModel<MC,RNG,S>::stats_type
        typedef typename MonteCarloModel<MC,RNG,S>::result_type
        virtual ~McSimulation() {}
        result_type value(Real tolerance,
                          Size maxSamples = QL_MAX_INTEGER,
                          Size minSamples = 1023) const;
        result_type valueWithSamples(Size samples) const;
        void calculate(Real requiredTolerance,
                       Size requiredSamples,
                       Size maxSamples) const;
        const stats_type& sampleAccumulator(void) const;
        McSimulation(bool antitheticVariate,
                     bool controlVariate);
        virtual shared_ptr<path_pricer_type> pathPricer() const = 0;
        virtual shared_ptr<path_generator_type> pathGenerator() const = 0;
        virtual TimeGrid timeGrid() const = 0;
        virtual shared_ptr<path_pricer_type>
        controlPathPricer() const;
        virtual shared_ptr<path_generator_type>
        controlPathGenerator() const;
        virtual result_type controlVariateValue() const;
        virtual shared_ptr<PricingEngine> controlPricingEngine() const;
        mutable shared_ptr<MonteCarloModel<MC,RNG,S> > mcModel_;
        bool antitheticVariate_, controlVariate_;

McSimulation follows the Template Method pattern. It asks the user to implement a few pieces of behavior, and in return it provides generic logic to instantiate a Monte Carlo model and run a simulation. Derived classes must define at least the pathPricer method, that returns the path pricer to be used; the pathGenerator method, that returns the path generator; and the timeGrid method, that returns the grid describing the nodes of the simulation. Other methods, returning the objects to be used for control variates, might or might not be McSimulation provides default implementations that return null values, so derived classes that don't want to use the control variate technique can just forget about it.

Based on such methods, McSimulation provides most of the behavior needed by an engine. Its constructor takes two boolean flags specifying whether it should use either antithetic or control variates; the second will only matter if the derived class implements the required methods.

The value method adds samples to the underlying model until the estimated error matches a required tolerance. (Of course, this needs an error estimate, so it won't work with low-discrepancy methods.)
    template <template <class> class MC, class RNG, class S>
    typename McSimulation<MC,RNG,S>::result_type
    McSimulation<MC,RNG,S>::value(Real tolerance,
                                  Size maxSamples,
                                  Size minSamples) const {
        Real error = mcModel_->sampleAccumulator().errorEstimate();
        while (error > tolerance) {
            QL_REQUIRE(sampleNumber < maxSamples, ...);
            Real order = (error*error)/(tolerance*tolerance);
            Size nextBatch =
            nextBatch = std::min(nextBatch, maxSamples-sampleNumber);
            sampleNumber += nextBatch;
            error = mcModel_->sampleAccumulator().errorEstimate();
        return mcModel_->sampleAccumulator().mean();
It looks at the current number of samples \( n \) and the current error \( \epsilon \), estimates the number of samples \( N \) that will be needed to reach the given tolerance \( \tau \) as \( N = n \times \epsilon^2/\tau^2 \) (since of course \( \epsilon \propto 1/\sqrt{n} \)), and adds a new batch of \( N-n \) samples that gets the total closer to the estimated number; then it assesses the error again, and repeats the process as needed.

The valueWithSamples method just adds a batch of samples so that their total number matches the required one; its implementation is not shown here because of space constraints, but is simple enough.

Finally, the calculate method runs a complete simulation.
    template <template <class> class MC, class RNG, class S>
    void McSimulation<MC,RNG,S>::calculate(Real requiredTolerance,
                                           Size requiredSamples,
                                           Size maxSamples) const {
        if (!controlVariate_) {
            mcModel_ = shared_ptr<MonteCarloModel<MC,RNG,S> >(
                new MonteCarloModel<MC,RNG,S>(
                           pathGenerator(), pathPricer(),
                           S(), antitheticVariate_));
        } else {
            ... // same as above, but passing the control-variate args, too.

        if (requiredTolerance != Null<Real>())
            this->value(requiredTolerance, maxSamples);
It takes as arguments either a required tolerance or a required number of samples (one of the arguments must be null, but not both) as well as a maximum number of samples; it instantiates a MonteCarloModel instance, with the controlVariate_ flag determining whether to pass the control-variate arguments to the constructor; and it calls either the value or the valueWithSamples method, depending on what goal was required.

In next post, I'll show an example of how to use the McSimulation class to build a pricing engine; but before that, let me point out a few ways in which it could be improved.

First of all, it currently implements logic to run a simulation based on two criteria (accuracy or total number of samples) but of course, more criteria could be added. For instance, one could run a simulation until a given clock time is elapsed; or again, one could add several batches of samples, look at the result after each one, and stop when convergence seems to be achieved. This suggests that the hard-coded value and valueWithSamples methods could be replaced by an instance of the Strategy pattern, and the switch between them in the calculate method by just a call to the stored strategy.

In turn, this would also remove the current ugliness in calculate: instead of passing the whole set of arguments for both calculations and giving most of them to a null value (like in the good old days, where languages could not overload methods) one would pass only the required arguments to the strategy object, and then the strategy object to calculate.

Finally, the presence of controlPricingEngine method is a bit of a smell. The implementation of McSimulation doesn't use it, so it's not strictly a part of the Template Method pattern and probably shouldn't be here. However, a couple of other classes (both inheriting from McSimulation, but otherwise unrelated) declare it and use it to implement the controlVariateValue method; therefore, leaving it here might not be the purest of designs but prevents some duplication.

Aside: synchronized walking

Both the MonteCarloModel and McSimulation class templates allow one to define an alternate path generator for control variates. Note, however, that this feature should be used with some caution. For this variance-reduction technique to work, the control-variate paths returned from the second generator must be correlated with the regular ones, which basically means that the two path generators must use two identical random-number generators: same kind, dimensionality, and seed (if any).

Unfortunately, this constraint rules out quite a few possibilities. For instance, if you're using a stochastic volatility model, such as the Heston model, you might be tempted to use the Black-Scholes model as control variate. No such luck: for \( n \) time steps, the Heston process requires \( 2n \) random numbers (\( n \) for the stock price and \( n \) for its volatility) while the Black-Scholes process just needs \( n \). This makes it impossible to keep the two corresponding path generators in sync.

In practice, you'll have a use case for the alternate path generator if you have a process with a number of parameters which is not analytically tractable in the generic case, but has a closed-form solution for your option value if some of the parameters are null. If you're so lucky, you can use a fully calibrated process to instantiate the main path generator, and another instance of the process with the null parameters to generate control-variate paths.

Liked this post? Share it: