Banner

Introduction to QuantLib Development with Luigi Ballabio, London, UK. October 19th to 21st, 2015.
Early-bird discount available until August 18th. Register now!

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:

Monday, June 15, 2015

Chapter 8, part 1 of n: the finite-differences framework.

Welcome back.

As I write, release candidates for QuantLib 1.6 have been out for a week or two (if you haven't heard about this already, you might want to subscribe to the QuantLib mailing list). If you have a few cycles to spare, please download it, give it a try and report any problems to the mailing list. The official 1.6 release should be out in a week or two. (Brace yourself for new developments afterwards. Just saying.)

This week, some initial content from chapter 8 of my book. Yes, I finally sat down to write it. The first couple subsections were pushed to the Leanpub edition a couple of weeks ago, with more coming during the summer (both there and 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.

The finite-differences framework

The last framework I'll be writing about was, in fact, the first to appear in the library. Of course, it was not a framework back then; it was just a bunch of code hard-coded inside a single option pricer.

Since then, the reusable parts of the code were reorganized and grew in complexity—probably to the point where it's hard to find one's way inside the architecture. A few years ago, a new version of the framework was contributed adding 2-dimensional finite-difference models (which were absent in the first version) and sporting a more modular architecture.

The two versions of the framework still coexist in the library. In this chapter, I'll describe both of them—while hoping for the new one to replace the first some day.

The old framework

As I mentioned in a previous post, our master plan to create a common framework for trees and finite-differences models never came to fruition. Instead, the finite-differences code stayed closer to its underlying mathematics and implemented concepts such as differential operators, boundary conditions and so on. I'll tackle one of them in each of the following subsections.

Differential operators

If I were to start at the very beginning, I'd describe the Array class here; but they behave as you can imagine, so you'll forgive me if I leave its description to a later post and just note that it's used here to discretize the function \( f(x) \) on a grid \( \{ x_0 \ldots x_{N-1} \} \) as the array \( \mathbf{f} \), with \( \mathbf{f}_i = f(x_i) \) for any \( i \).

Onwards to operators. Now, I'm probably preaching to the choir, so I'll keep the math to a minimum to avoid telling you lots of stuff you know already. (And if you don't, there's no way I can tell you all you need. You can refer to any of the several books on finite differences; for instance, Duffy, 2006 [1]).

In short: a differential operator transforms a function \( f(x) \) in one of its derivatives, say, \( f^{\prime}(x) \). The discretized version transforms \( \mathbf{f} \) into \( \mathbf{f^{\prime}} \), and since differentiation is linear, it can be written as a matrix. In general, the operator doesn't give the exact discretization of the derivative but just an approximation; that is, \( \mathbf{f^{\prime}}_i = f^{\prime}(x_i) + \epsilon_i \) where the error terms can be reduced by decreasing the spacing of the grid.

QuantLib provides differential operators for the first derivative \( \partial/\partial x \) and the second derivative \( \partial^2/\partial x^2 \). They are called \( D_0 \) and \( D_+D_- \), respectively, and are defined as:
$$
\frac{\partial f}{\partial x}(x_i) \approx D_0 \mathbf{f}_i
= \frac{\mathbf{f}_{i+1}-\mathbf{f}_{i-1}}{2h} \qquad
\frac{\partial^2 f}{\partial x^2}(x_i) \approx D_+D_- \mathbf{f}_i
= \frac{\mathbf{f}_{i+1}-2\mathbf{f}_i+\mathbf{f}_{i-1}}{h^2}
$$
where \( h = x_i - x_{i-1} \) (we're assuming that the grid is evenly spaced). Taylor expansion shows that the error is \( o(h^2) \) for both of them.

As you can see, the value of the derivative at any given index \( i \) (except the first and last, that I'll ignore now) only depends from the values of the function at the same index, the one that precedes it, and the one that follows. This makes \( D_0 \) and friends tridiagonal: the structure of such operators can be sketched as
$$
\left(
\begin{array}{ccccccc}
d_o & u_0 & & & & & \\
l_0 & d_1 & u_1 & & & & \\
& l_1 & d_2 & u_2 & & & \\
& & \ddots & \ddots & \ddots & & \\
& & & l_{n-4} & d_{n-3} & u_{n-3} & \\
& & & & l_{n-3} & d_{n-2} & u_{n-2} \\
& & & & & l_{n-2} & d_{n-1}
\end{array}
\right)
$$

They have a number of desirable properties, including the fact that a linear combination of tridiagonal operators (such as the one found, say, in the Black-Scholes-Merton formula) is still tridiagonal; thus, instead of using a generic Matrix class (which does exist, but is not used here), we made them instances of a specific class called TridiagonalOperator and shown in the following listing.
    class TridiagonalOperator {
        Size n_;
        Array diagonal_, lowerDiagonal_, upperDiagonal_;

      public:
        typedef Array array_type;

        explicit TridiagonalOperator(Size size = 0);
        TridiagonalOperator(const Array& low,
                            const Array& mid,
                            const Array& high);

        static Disposable<TridiagonalOperator>
        identity(Size size);

        Size size() const;
        const Array& lowerDiagonal() const;
        const Array& diagonal() const;
        const Array& upperDiagonal() const;

        void setFirstRow(Real, Real);
        void setMidRow(Size, Real, Real, Real);
        void setMidRows(Real, Real, Real);
        void setLastRow(Real, Real);

        Disposable<Array> applyTo(const Array& v) const;
        Disposable<Array> solveFor(const Array& rhs) const;
        void solveFor(const Array& rhs, Array& result) const;

      private:
        friend Disposable<TridiagonalOperator>
        operator+(const TridiagonalOperator&,
                  const TridiagonalOperator&);
        friend Disposable<TridiagonalOperator>
        operator*(Real, const TridiagonalOperator&);
        // other operators
    };
As you can see from the private members, the representation of the operator reflects its tridiagonal structure: instead of storing all the mostly null elements of the matrix, we just store three arrays corresponding to the diagonal, lower diagonal, and upper diagonal; that is, the \( d_i \), \( l_i \) and \( u_i \) of the earlier sketch.

The constructors and the first bunch of methods deal with building and inspecting the operator. The first constructor sets the size of the matrix (that is, of the three underlying arrays) and trusts the values to be filled later, by means of the provided setters: referring to the figure again, the setFirstRow method sets \( d_0 \) and \( u_0 \), the setLastRow method sets \( l_{n-2} \) and \( d_{n-1} \), setMidRow sets \( l_{i-1} \), \( d_i \) and \( u_i \) for a given \( i \), and the convenience method setMidRows sets them for all \(i \). The second constructor sets all three arrays in one fell swoop; and a factory method, identity, works as a third constructor by building and returning an operator with \( 1 \) on the diagonal and \( 0 \) elsewhere. The inspectors do the obvious thing and return each the corresponding data member. Oftentimes, the returned values are wrapped in a Disposable class template, in a possibly misguided attempt to avoid a few copies. Again, details will be in a future post.

The applyTo and solveFor methods are the meat of the interface. Given a tridiagonal operator L and an array u, L.applyTo(u) returns the result \( L u \) of applying \( L \) to \( u \), while L.solveFor(u) returns the array \( v \) such that \( u = L v \) (The overloaded call L.solveFor(u,v) does the same while avoiding an allocation). Both operations will be used later, and for a tridiagonal operator they both have the nice property to be \( O(N) \) in time, that is, to only take time proportional to the size \( N \) of the operator. The algorithms were taken from the classic Numerical Recipes in C [2]. Not the code, of course: besides not being free, the original code had a number of endearing old conventions (such as numbering arrays from 1, like in Downton Abbey) that wouldn't have matched the rest of the library.

Finally, the class defines a few operators. They make it possible to write numerical expressions such as 2*L1 + L2 that might be used to compose tridiagonal operators. To this effect, the library also provides a few building blocks such as the \( D_0 \) and \( D_+D_- \) operators I mentioned earlier; for instance, given the Black-Scholes equation, written in operator form as
$$
\frac{\partial f}{\partial t} =
\left[
- \left(r - q - \frac{\sigma^2}{2} \right) \frac{\partial}{\partial x}
- \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2}
+ r I
\right] f \equiv L_{BS} f
$$
one could define the corresponding \( L_{BS} \) operator as
    TridiagonalOperator L_BS = -(r-q-sigma*sigma/2)*DZero(N)
                               -(sigma*sigma/2)*DPlusDMinus(N)
                               -r*TridiagonalOperator::identity(N);
In this case, though, we didn't eat our own dog food. The library does define a Black-Scholes operator, but not through composition as above.

A final note: while they're probably the most performant, tridiagonal operators are not the only ones that could be used in a finite-difference scheme. However, there's no base class for such operators. The interface that they should implement (that is, applyTo and solveFor is only defined implicitly, through its use in the template classes that we'll see in the next subsection. In hindsight, defining such a base class wouldn't have degraded performance; the cost of a virtual call to applyTo or solveFor would have been negligible, compared to the cost of the actual calculation. But we were eagerly reading about template numerical techniques at that time; and like the two guys in the Wondermark strip, we decided to jingle all the way.

Bibliography

[1] D. J. Duffy, Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach. John Wiley and Sons, 2006.
[2] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd edition. Cambridge University Press, 1992.

Liked this post? Share it: