Banner

Introduction to QuantLib Development with Luigi Ballabio, London, UK. October 19th to 21st, 2015.
Group discounts available. Register now!

Monday, August 24, 2015

Chapter 8, part 5 of n: finite-difference models

Welcome back, everybody.

This week, the next part of the series on the finite-difference framework.

But first, a small personal milestone: while I was in vacation, the 100th copy of Implementing QuantLib was sold on Leanpub. So, to the 100th reader who bought it between last Friday and Saturday, but also to all the others: thank you.

A couple of weeks ago we also had... well, actually, I don't know if we had some kind of event or not. One morning I found in my inbox a few hundreds notifications that email addresses were removed from the QuantLib mailing lists. This usually happens when an address bounces a few times; the mailing-list software then flags it as obsolete and removes it. However, we never had that many at once. I'm not sure if this is some kind of glitch after the recent Sourceforge outages, or the result of their changing some kind of threshold. In any case, if you haven't received emails from the list in a while, please check that you're still subscribed. You can do so from http://quantlib.org/mailinglists.shtml: follow the "Subscribe/Unsubscribe/Preferences" link for the lists you're subscribed to.

In the "back to business" department: there's a few things happening in October. If you're going to the WBS Fixed Income Conference in Paris, you may find me there. On Wednesday 7trh, I'll be in the bench during Alexander Sokol's workshop on how to add AAD to QuantLib (I'm listed as a presenter, but I'll be just there ready to intervene if bits of particular QuantLib knowledge are needed). On the morning of Friday 9th, instead, I'll have a slot to present QuantLib and its next developments in fixed income.

Also in October, I'll be running my next Introduction to QuantLib Development course (in London, from Monday 19th to Wednesday 21st). Details and a registration form are at http://bit.ly/quantlib2015; you can also check out the hangout I did with Jacob Bettany to talk about it.

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

The FiniteDifferenceModel class

The FiniteDifferenceModel class, shown in the next listing, uses an evolution scheme to roll an asset back from a later time \( t_1 \) to an earlier time \( t_2 \), taking into account any conditions we add to the model and stopping at any time we declare as mandatory (because something is supposed to happen at that point).

A word of warning: as I write this, I'm looking at this class after quite a few years and I realize that some things might have been done differently—not necessarily better, mind you, but there might be alternatives. After all this time, I can only guess at the reasons for such choices. While I describe the class, I'll point out the alternatives and write my guesses.
    template<class Evolver>
    class FiniteDifferenceModel {
      public:
        typedef typename Evolver::traits traits;
        typedef typename traits::operator_type operator_type;
        typedef typename traits::array_type array_type;
        typedef typename traits::bc_set bc_set;
        typedef typename traits::condition_type condition_type;

        FiniteDifferenceModel(
            const Evolver& evolver,
            const std::vector<Time>& stoppingTimes =
                                          std::vector<Time>());

        FiniteDifferenceModel(
            const operator_type& L,
            const bc_set& bcs,
            const std::vector<Time>& stoppingTimes =
                                          std::vector<Time>());
        const Evolver& evolver() const{ return evolver_; }
        void rollback(array_type& a,
                      Time from,
                      Time to,
                      Size steps) {
             rollbackImpl(a,from,to,steps,
                          (const condition_type*) 0);
        }
        void rollback(array_type& a,
                      Time from,
                      Time to,
                      Size steps,
                      const condition_type& condition) {
            rollbackImpl(a,from,to,steps,&condition);
        }
      private:
        void rollbackImpl(array_type& a,
                          Time from,
                          Time to,
                          Size steps,
                          const condition_type* condition);
        Evolver evolver_;
        std::vector<Time> stoppingTimes_;
    };
I'll just mention quickly that the class takes the evolution scheme as a template argument (as in the code, I might call it evolver for brevity from now on) and extracts some types from it via traits; that's a technique you saw a few times already in earlier chapters.

So, first of all: the constructors. As you see, there's two of them. One, as you'd expect, takes an instance of the evolution scheme and a vector of stopping times (and, not surprisingly, stores copies of them). The second, instead, takes the raw ingredients and builds the evolver.

The reason for the second one might have been convenience, that is, to avoid writing the slighly redundant
    FiniteDifferenceModel<SomeScheme>(SomeScheme(L, bcs), ts);
but I'm not sure that the convenience offsets having more code to maintain (which, by the way, only works for some evolvers and not for others; for instance, the base MixedScheme class doesn't have a constructor taking two arguments). This constructor might make more sense when using one of the typedefs provided by the library, such as StandardFiniteDifferenceModel, which hides the particular scheme used; but I'd argue that you should know what scheme you're using anyway, since they have different requirements in terms of grid spacing and time steps. (By the same argument, the StandardFiniteDifferenceModel typedef would be a suggestion of what scheme to use instead of a way to hide it.)

The main logic of the class is coded in the rollback and rollbackImpl methods. The rollback method is the public interface and has two overloaded signatures: one takes the array of values to roll back, the initial and final times, and the number of steps to use, while the other also takes a StepCondition instance to use at each step. Both forward the call to the rollbackImpl method, which takes a raw pointer to step condition as its last argument; the pointer is null in the first case and contains the address of the condition in the second.

Now, this is somewhat unusual in the library: in most other cases, we'd just have a single public interface taking a shared_ptr and defining a null default value. I'm guessing that this arrangement was made to allow client code to simply create the step condition on the stack and pass it to the method, instead of going through the motions of creating a shared pointer (which, in fact, only makes sense if it's going to be shared; not passed, used and forgotten). Oh, and by the way, the use of a raw pointer in the private interface doesn't contradict the ubiquitous use of smart pointers everywhere in the library. We're not allocating a raw pointer here; we're just passing the address of an object living on the stack (or inside a smart pointer, for that matters) and which can manage its lifetime very well, thank you.

One last remark on the interface: we're treating boundary conditions and step conditions differently, since the former are passed to the constructor of the model and the latter are passed to the rollback method. (Another difference is that we're using explicit sets of boundary conditions but just one step condition, which can be a composite of several ones if needed. I have no idea why.) One reason could have been that this allows more easily to use different step conditions in different calls to rollback, that is, at different times; but that could be done inside the step condition itself, as it is passed the current time. The real generalization would have been to pass the boundary conditions to rollback, too, as it's currently not possible to change them at a given time in the calculation. This would have meant changing the interfaces of the evolvers, too: the boundary conditions should have been passed to their step method, instead of their constructors.

And finally, the implementation.
    void FiniteDifferenceModel::rollbackImpl(
                               array_type& a,
                               Time from,
                               Time to,
                               Size steps,
                               const condition_type* condition) {
        Time dt = (from-to)/steps, t = from;
        evolver_.setStep(dt);
        if (!stoppingTimes_.empty() && stoppingTimes_.back()==from) {
            if (condition)
                condition->applyTo(a,from);
        }
        for (Size i=0; i<steps; ++i, t -= dt) {
            Time now = t, next = t-dt;
            if (std::fabs(to-next) < std::sqrt(QL_EPSILON))
                next = to;
            bool hit = false;
            for (Integer j = stoppingTimes_.size()-1; j >= 0 ; --j) {
                if (next <= stoppingTimes_[j]
                         && stoppingTimes_[j] < now) {
                    hit = true;
                    evolver_.setStep(now-stoppingTimes_[j]);
                    evolver_.step(a,now);
                    if (condition)
                        condition->applyTo(a,stoppingTimes_[j]);
                    now = stoppingTimes_[j];
                }
            }
            if (hit) {
                if (now > next) {
                    evolver_.setStep(now - next);
                    evolver_.step(a,now);
                    if (condition)
                        condition->applyTo(a,next);
                }
                evolver_.setStep(dt);
            } else {
                evolver_.step(a,now);
                if (condition)
                    condition->applyTo(a, next);
            }
        }
    }
We're going back from a later time from to an earlier time to in a given number of steps. This gives us a length dt for each step, which is set as the time step for the evolver. If from itself is a stopping time, we also apply the step condition (if one was passed; this goes for all further applications). Then we start going back step by step, each time going from the current time t to the next time, t-dt; we'll also take care that the last time is actually to and possibly snap it back to the correct value, because floating-point calculations might have got us almost there but not quite.

If there were no stopping times to hit, the implementation would be quite simple; in fact, the whole body of the loop would be the three lines inside the very last else clause:
    evolver_.step(a,now);
    if (condition)
        condition->applyTo(a, next);
that is, we'd tell the evolver to roll back the values one step from now, and we'd possibly apply the condition at the target time next. This is what actually happens when there are no stopping times between now and next, so that the control variable hit remains false.

In case there are one or more stopping times in the middle of the step, things get a bit more complicated. Have a look at the upper half of the following figure, where I sketched a series of regular steps with two stopping times, \( s_1 \) and \( s_2 \), falling in the middle of two of them.


Let's say we have just made the step from \( t_7 \) to \( t_6 \) as described in the previous paragraph, and should now go back to \( t_5 \). However, the code will detect that it should stop at \( s_2 \), which is between now and next; therefore, it will set the control variable hit to true, it will set the step of the evolver to the smaller step \( t_6-s_2 \) that will bring it to the stopping time, and will perform the step (again, applying the condition afterwards). Another, similar change of step would happen if there were two or more stopping times in a single step; the code would then step from one to the other. Finally, the code will enter the if (hit) clause, which get things back to normal: it performs the remaining step \( s_2-t_5 \) and then resets the step of the evolver to the default value, so that it's ready for the next regular step to \( t_4 \).

As a final remark, let's go back to trees for a minute. Talking about the tree framework, I mentioned in passing that the time grid of a lattice is built to include the mandatory stopping times. In this case, the grid would have been the one shown in the lower half of the previous figure. As you see, the stopping times \( s_1 \) and \( s_2 \) were made part of the grid by using slightly smaller steps between to and \( s_1 \) and slightly larger steps between \( s_2 \) and from. (In this case, the change results in the same number of steps as the finite-difference grid. This is not always guaranteed; the tree grid might end up having a couple of steps more or less than the requested number.) I'm not aware of drawbacks in using either of the two methods, but experts might disagree so I'm ready to stand corrected. In any case, should you want to replicate the tree grid in a finite-difference model, you can do so by making several calls to rollback and going explicitly from stopping time to stopping time.

Liked this post? Share it:

Monday, July 27, 2015

Chapter 8, part 4 of n: step conditions

Welcome back.

The content for this post is the fourth part of an ongoing series on the QuantLib finite-difference framework.

But first: this week the final version of Visual Studio 2015 was released and, whaddayaknow, QuantLib 1.6 was not quite prepared for it. We had tried to be forward-looking and had prepared our version-detecting code for VC13, but alas, Visual Studio 2015 turned out to be VC14 instead. The library compiles and links just fine, but has the wrong name and the auto-link pragmas get confused. This week I might be releasing a 1.6.1 version that will fix the issue, God willing and the creek don't rise—since the other bit of news this week was that Sourceforge had a major outage (you might have noticed that the QuantLib web site and the mailing lists have been down for a few days). Now it's mostly up again, but as of today it's still not possible to make new releases. Keep your fingers crossed.

Oh, and my Quants Hub workshop is still a steal at £99, but from what I'm told this will only last for a few more days, until the end of July. Get it cheap while you can.

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

Step conditions

At any step (or in the middle of a step; we'll see later how we account for this) something might happen that changes the price of the asset; for instance, a right can be exercised or a coupon could be paid.

As opposed to boundary conditions, we called these step conditions. The name might not be the best; I can easily see the model as enforcing the "condition" that the asset price is the maximum between its continuation and the intrinsic value, but a coupon payment is harder to categorize in this way. In any case, the base class for all such conditions is shown in the listing below.
    template <class array_type>
    class StepCondition {
      public:
        virtual ~StepCondition() {}
        virtual void applyTo(array_type& a, Time t) const = 0;
    };
The applyTo method, which modifies the array of asset values in place, must also take care of checking that the condition applies at the given time; so, for instance, an American exercise will change the asset values regardless of the time, whereas a dividend payment will first check the passed time and bail out if it doesn't match the known payment time.

As you see, the base class is quite simple; but derived classes over-generalized and went South. In the next listing, you can see one such class, from which even simple conditions (like the one for American exercise) are derived. We could have implemented a simple applyTo method comparing the option values and the intrinsic values. Instead, we started adding constructors over the years so that instances could be built from either an existing array of intrinsic values, or a payoff, or some data that could describe the payoff. Of course these alternatives needed to be stored into different data members, so we added type erasure to the mix, we defined an interface to retrieve intrinsic values, and we gave it an array-based and a payoff-based implementation.
    template 
    class CurveDependentStepCondition
        : public StepCondition {
      public:
        void applyTo(Array &a, Time) const {
            for (Size i = 0; i < a.size(); i++) {
                a[i] = applyToValue(a[i], getValue(a,i));
            }
        }
      protected:
        CurveDependentStepCondition(Option::Type type,
                                    Real strike);
        CurveDependentStepCondition(const Payoff *p);
        CurveDependentStepCondition(const array_type & a);

        class CurveWrapper {
          public:
            Real getValue(const array_type &a,
                          Size index) const = 0;
        };
        boost::shared_ptr curveItem_;

        Real getValue(const array_type &a, Size index) const {
            return curveItem_->getValue(a, index);
        }

        virtual Real applyToValue(Real, Real) const = 0;

        class ArrayWrapper : public CurveWrapper {
            array_type value_;
          public:
            ...
        };

        class PayoffWrapper : public CurveWrapper {
            boost::shared_ptr payoff_;
          public:
            ...
        };
    };
Ironically, we only ever used half of this. Existing engines use the array-based constructor and implementation, and the payoff-based inner class is not used anywhere in the library, Which is just as well; because, as I reviewed the code to write this chapter, I found out that its implementation is wrong. So all in all, we could have written a much simpler implementation of an American exercise condition that would just take an array of intrinsic values and implement applyTo accordingly. If you ever implement a step condition, I suggest you do the same.

A final note; the step condition is not applied inside the step method of the evolution scheme, as we'll see shortly. What this means is that there's no guarantee that the step condition won't break the boundary condition, or even the other way around. In practice, though, applying a step condition for a model that makes any sense will naturally enforce the chosen boundary condition.

Next time, with all the basic pieces in place, we'll start to put them together.

Aside: look, Ma, no hands

The faulty implementation of PayoffWrapper was added in February 2003, which probably makes this the longest-lived bug in the library (apart from those we still haven't found, of course). If I hadn't been looking at the code to write this chapter, the critter could have got old enough to drive.

Now, bugs happen; nothing shocking about this. What irks me, though, is that the problem would have been caught easily by a couple of unit tests. I and others usually try and review what goes in the library, and the code that gets contributed has often been used in the field and purged of obvious errors; but still I'm somewhat worried about the lack of coverage of old code in the test suite. If you have a coverage tool and some time on your hands, please let me hear from you. I'd love to see some analysis on this.

Liked this post? Share it:

Monday, July 20, 2015

MoneyScience Hangout available on YouTube

Hello everybody.

A short bit of news: that Google Hangout I told you I would do with MoneyScience's Jacob Bettany? We did it, and it's now online on YouTube. Here it is. More details on my next course (which will be in London on October 19th to 21st) are at http://bit.ly/quantlib2015.



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

Liked this post? Share it:

Monday, July 13, 2015

Chapter 8, part 3 of n: boundary conditions

Hello everybody.

Just a reminder before this week's content (that is, another installment in the series on finite-difference methods): my Quants Hub workshop, A Look at QuantLib Usage and Development, is going to be £99 until the end of July. Grab it cheap while you can.

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.

Boundary conditions

As I said, there's more going on in the step method. To begin with, we'll have to enforce boundary conditions at either end of the grid; for instance, we might want to ensure that an option has price close to 0 at the grid boundary which is deep out of the money and derivative close to 1 at the boundary which is deep in the money. Such conditions are modeled by the BoundaryCondition class, shown in the listing below.
    template <class Operator>
    class BoundaryCondition {
      public:
        typedef Operator operator_type;
        typedef typename Operator::array_type array_type;

        enum Side { None, Upper, Lower };

        virtual ~BoundaryCondition() {}

        virtual void applyBeforeApplying(operator_type&) const = 0;
        virtual void applyAfterApplying(array_type&) const = 0;
        virtual void applyBeforeSolving(operator_type&,
                                        array_type& rhs) const = 0;
        virtual void applyAfterSolving(array_type&) const = 0;
    };
As you see, the class uses runtime polymorphism. This is because we wanted to store collections of them, possibly of different types (you might remember the bc_set typedef from the listing of the MixedScheme class), and most suitable containers are homogeneous and thus require a common base class (std::tuple wasn't even a gleam in the standard committee's eye, and handling it requires template metaprogramming, which I'd try to avoid unless necessary—even in this time and age) We could have used std::pair if we had decided to limit ourselves to the 1-D case; we didn't, though, and in hindsight I think it avoided adding more complexity to the framework.

The interface of the class is, well, peculiar. Glossing over the Side enumeration (which is meant to specify a grid side in 1-D models, and thus too specific), methods like applyBeforeApplying are plausible contenders for the worst name in the library. The idea here is that the boundary condition can be enforced, or "applied", in different ways. One can wait for the operator's applyTo or solveFor to execute, and then modify the result; this would be implemented by applyAfterApplying and applyAfterSolving, respectively. Another possibility is to set up the operator and the input array beforehand, so that the result of the operation will satisfy the condition; this is done in applyBeforeApplying and applyBeforeSolving. (For some reason, the former only takes the operator. It might have been an oversight, due to the fact that the existing conditions didn't require to modify the array.)

As an example, look at the DirichletBC class in the listing below. It implements a simple Dirichlet boundary condition, i.e., one in which the function value at a given end of the grid must equal a given constant value.
    class DirichletBC
        : public BoundaryCondition<TridiagonalOperator> {
      public:
        DirichletBC(Real value, Side side);

        void applyBeforeApplying(TridiagonalOperator& L) const {
            switch (side_) {
              case Lower:
                L.setFirstRow(1.0,0.0);
                break;
              case Upper:
                L.setLastRow(0.0,1.0);
                break;
            }
        }
        void applyAfterApplying(Array& u) const {
            switch (side_) {
              case Lower:
                u[0] = value_;
                break;
              case Upper:
                u[u.size()-1] = value_;
                break;
            }
        }
        void applyBeforeSolving(TridiagonalOperator& L,
                                Array& rhs) const {
            switch (side_) {
              case Lower:
                L.setFirstRow(1.0,0.0);
                rhs[0] = value_;
                break;
              case Upper:
                L.setLastRow(0.0,1.0);
                rhs[rhs.size()-1] = value_;
                break;
            }
        }
        void applyAfterSolving(Array&) const {}
    };
As you can see, it is no longer a class template; when implemented, boundary conditions are specialized for a given operator—not surprisingly, because they have to know how to access and modify it. In this case, we're using the tridiagonal operator I described in an earlier subsection.

The constructor is simple enough: it takes the constant value and the side of the grid to which it must be applied, and stores them.

As I said, the applyBeforeApplying method must set up the operator L so that the result of L.apply(u), or \( \mathbf{u^{\prime}} = L \cdot \mathbf{u} \), satisfies the boundary condition. To do this, it sets the first (last) row to the corresponding row of the identity matrix; this ensures that the first (last) element of \( \mathbf{u^{\prime}} \) equal the corresponding element of \( \mathbf{u} \). Given that the array satisfied the boundary condition at the previous step, it keeps satisfying it at this step. Relying on this is a bit unsafe (for instance, it would break if the value changed at different times, or if any other event modified the array values) but it's the best one can do without accessing the input array. The applyAfterApplying method is simpler and safer: it just sets the value of the output array directly.

The applyBeforeSolving method works as its sibling applyBeforeApplying, but it can also access the input array, so it can ensure that it contains the correct value. Finally, and surprisingly enough, as of this writing the applyAfterSolving method is empty. When we wrote this class, we probably relied on the fact that applyBeforeSolving would also be called, which doesn't strike me as a very good idea in hindsight: we're not sure of what methods any given evolution scheme may call (this also makes applyBeforeApplying less safe than it seems).

The reason this worked is the way the MixedScheme::step method is implemented, as seen in the next listing. At each step, all boundary conditions are enforced first before and after applying the explicit part of the scheme, and then before and after solving for the implicit part. In the case of the Dirichlet condition, this causes applyAfterApplying to fix any problem that applyBeforeApplying might have introduced and also ensures that applyBeforeSolving does the work even if applyAfterSolving doesn't.

If a boundary condition class implemented all of its methods correctly, this would be redundant; enforcing the condition just after the operations would be enough. However, we probably assumed that some conditions couldn't (or just didn't) implement all of them, and went for the safe option.
    template <class Operator>
    void MixedScheme<Operator>::step(array_type& a, Time t) {
        if (theta_ != 1.0) { // there is an explicit part
            for (Size i=0; i<bcs_.size(); i++)
                bcs_[i]->applyBeforeApplying(explicitPart_);
            a = explicitPart_.applyTo(a);
            for (Size i=0; i<bcs_.size(); i++)
                bcs_[i]->applyAfterApplying(a);
        }
        if (theta_ != 0.0) { // there is an implicit part
            for (Size i=0; i<bcs_.size(); i++)
                bcs_[i]->applyBeforeSolving(implicitPart_, a);
            implicitPart_.solveFor(a, a);
            for (Size i=0; i<bcs_.size(); i++)
                bcs_[i]->applyAfterSolving(a);
        }
    }
If we ever set out to revamp this part of the library (not that likely, given the new framework) this is something we should get a look at. The problem to solve, unfortunately, is not an easy one: to call only the methods that are required by the scheme and supported by the boundary condition, we'd need to know the specific type of both.

Inverting the dependency and passing the scheme to the boundary condition, instead of the other way around, wouldn't work: the boundary condition wouldn't know if the scheme is calling applyTo or solveFor, and wouldn't be able to act inside the scheme's step method (as in MixedScheme::step, where the boundary conditions are enforced between the implicit and explicit parts of the step.

A mediator or some kind of visitor might work, but the complexity of the code would increase (especially if the number of conditions grow). Had they been available in C++, multimethods might have been the best choice; but even those might not be usable for multiple conditions. All in all, the current workaround of calling all boundary condition method might be redundant but at least works correctly. The only alternative I can think of would be to require all conditions to be enforced after the operation, remove the applyBefore... methods and be done with it.

Enough with boundary conditions for now. With the pieces we have so far, we can take a payoff at \( t=T \) and evolve it back step by step to \( t=0 \)—if nothing happens in between, that is. But that's not often the case, which is my cue for the next post.

Liked this post? Share it:

Monday, July 6, 2015

PSA: Google Hangout

Hello everybody.

A very short post for a public service announcement. The folks at MoneyScience have organized a Google Hangout with me on next Monday, July 13th. Details are at this link.

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

Liked this post? Share it:

Monday, June 29, 2015

A quick look at the QuantLib 1.6 release

Hello everybody.

Unlike in Westeros, summer is coming.  My posting schedule will be more relaxed, as usual, so don't hold your breath waiting for the next posts (hint: subscribe instead, and you'll be notified whenever I post).

A bit of news: the guys at Quants Hub are happy of the success of their new pricing structure, and they're extending it to July 31st. This means that, for the whole next month, you can still get my video workshop for £99.

However, the bigger piece of news is that last Tuesday I've released QuantLib 1.6; you can download it at this link if you haven't already. Thanks to all those who contributed. In this post, I'll share some information on the release (like I did last time).

I'm happy to report that we managed to decrease the distance between releases. This one comes out four and a half months after QuantLib 1.5, which is much better than the one-year hiatus between the 1.4 and 1.5 releases. However, this doesn't mean that the 1.6 release was a small one: it contains 65 issues and pull requests (follow the link to see their list and examine the changes brought by any one of them). A couple of quick git commands also shows a grand total of 382 commits by 16 contributors (git outputs 18 because it includes a couple of duplicates it doesn't know about).


As usual, the list doesn't include contributors that reported bugs, suggested changes, or even contributed patches but didn't actually author the relevant commits. Their names are in the list of changes for the 1.6 release; apologies if I missed anyone.

As I write, there's 27 open pull requests already in the queue for next release. Stay tuned.

Again, thanks to all those who contributed to the 1.6 release!

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

Liked this post? Share it:

Monday, June 22, 2015

Chapter 8, part 2 of n: evolution schemes

Welcome back.

This week, the second part of the series on the finite-differences framework that started in the last post. And also this week, I'll be releasing QuantLib 1.6. Stay tuned.

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.

Evolution schemes

In a partial differential equation \( {\partial f}/{\partial t} = L f \), the operators in the previous subsection provide a discretization of the derivatives on the right-hand side. Evolution schemes discretize the time derivative on the left-hand side instead.

As you know, finite-difference methods in finance start from a known state \( \mathbf{f}(T) \) at maturity (the payoff of the derivative) and evolve backwards to today's date. At each step, we need to evaluate \( \mathbf{f}(t) \) based on \( \mathbf{f}(t+\Delta t) \).

The simplest way is to approximate the equation above as
$$
\frac{\mathbf{f}(t+\Delta t) - \mathbf{f}(t)}{\Delta t} =
L \cdot \mathbf{f}(t+\Delta t)
$$
which simplifies to \( \mathbf{f}(t) = \mathbf{f}(t+\Delta t) - \Delta t \cdot L \cdot \mathbf{f}(t+\Delta t) \), or \( \mathbf{f}(t) = (I - \Delta t \cdot L) \cdot \mathbf{f}(t+\Delta t) \). It's a simple matrix multiplication, and updating the function array translates into
    a = (I - dt*L).applyTo(a);
(in actual code, of course, the I + dt*L operator would be precalculated and stored across time steps).

A slightly less simple way is to write
$$
\frac{\mathbf{f}(t+\Delta t) - \mathbf{f}(t)}{\Delta t} =
L \cdot \mathbf{f}(t)
$$
which results into \( (I + \Delta t \cdot L) \cdot \mathbf{f}(t) = \mathbf{f}(t + \Delta t) \). This is a more complex step, since it requires solving a linear system (or inverting the operator, which is the same); but given the interface of our operators, the update translates equally simply into
    a = (I + dt*L).solveFor(a);
Now, this is where using tridiagonal operators pays off. Let me oversimplify, so I can keep a long story short: schemes of the first kind (called explicit schemes) are simpler but less stable. Implicit schemes, instead (those of the second kind), are more stable and can be used with larger time steps. If the operator were a generic one, the price for an implicit step would be \( O(N^3) \) in the size of the grid, making it unpractical when compared to the \( O(N^2) \) cost for an explicit step. For tridiagonal operator, though, both methods cost \( O(N) \), which allows us to reap the benefits of explicit steps with no additional price.

Going back to the library: the two schemes above (called explicit Euler and implicit Euler scheme, respectively) are both available. They are two particular specializations of a generic MixedScheme class template, shown in the listing that follows, which models the discretization scheme
$$
\frac{\mathbf{f}(t+\Delta t) - \mathbf{f}(t)}{\Delta t} =
L \cdot \left[(1-\theta) \cdot \mathbf{f}(t+\Delta t)
+ \theta \cdot \mathbf{f}(t) \right]
$$
where an implicit and explicit step are mixed. As you can see from the formula, \( \theta = 0 \) gives the explicit Euler scheme and \( \theta = 1 \) gives implicit Euler; any other value in the middle can be used (for instance, \( 1/2 \) gives the Crank-Nicolson scheme).
    template <class Operator>
    class MixedScheme  {
      public:
        typedef OperatorTraits<Operator> traits;
        typedef typename traits::operator_type operator_type;
        typedef typename traits::array_type array_type;
        typedef typename traits::bc_set bc_set;
        typedef typename traits::condition_type condition_type;

        MixedScheme(const operator_type& L,
                    Real theta,
                    const bc_set& bcs)
        : L_(L), I_(operator_type::identity(L.size())),
          dt_(0.0), theta_(theta) , bcs_(bcs) {}

        void setStep(Time dt) {
            dt_ = dt;
            if (theta_!=1.0) // there is an explicit part
                explicitPart_ = I_-((1.0-theta_) * dt_)*L_;
            if (theta_!=0.0) // there is an implicit part
                implicitPart_ = I_+(theta_ * dt_)*L_;
        }
        void step(array_type& a, Time t) {
            if (theta_!=1.0) {
                a = explicitPart_.applyTo(a);
            if (theta_!=0.0) {
                implicitPart_.solveFor(a, a);
        }
      protected:
        operator_type L_, I_, explicitPart_, implicitPart_;
        Time dt_;
        Real theta_;
        bc_set bcs_;
    };

    template <class Operator>
    class ExplicitEuler : public MixedScheme<Operator> {
      public:
        // typedefs, not shown
        ExplicitEuler(const operator_type& L,
                      const bc_set> >& bcs)
        : MixedScheme<Operator>(L, 0.0, bcs) {}
    };
As I mentioned before, schemes take the type of their differential operators as a template argument. The MixedScheme class does that as well, and extracts a number of types from the operator class by means of the OperatorTraits template (that I'll describe later on).

The constructor of the MixedScheme class takes the differential operator \( L \), the value of \( \theta \), and a set of boundary conditions (more on that later) and stores them. Another method, setStep, stores the \( \Delta t \) to be used as time step and precomputes the operators \( I - (1-\theta) \cdot \Delta t \cdot L \) and \( I + \theta \cdot \Delta t \cdot L \) on which the applyTo and solveFor methods will be called; the two special cases \( \theta = 0 \) and \( \theta = 1 \) are checked in order to avoid unnecessary computations. This is done in a separate method in case we needed to change the time step in the middle of a simulation; by being able to reset the time step, we can reuse the same scheme instance.

Finally, the step method performs the actual evolution of the function array. The implementation shown in the listing is just a basic sketch that I'll be improving in the next subsections; here, it just performs the explicit and implicit part of the scheme in this order (again, avoiding unnecessary operations).

The few schemes provided by the library (such as the ExplicitEuler, shown in the listing) are all specializations of the MixedScheme class (I'm just talking about the old framework here; the new framework has different ones). Again, there's no base class for an evolution scheme; other parts of the framework will take them as a template argument and expect them to provide the setStep and step methods.

Liked this post? Share it: