Welcome back.

This week’s content is the fifth part of the series on the QuantLib tree framework which started in this post.

Subscribe to my Substack to receive my posts in your inbox, or follow me on Twitter or LinkedIn if you want to be notified of new posts, or subscribe via RSS if you’re the tech type: the buttons for all that are in the footer. Also, I’m available for training, both online and (when possible) on-site: visit my Training page for more information.

The TreeLattice class template

After the last couple of posts on trees, we’re now back to lattices. The TreeLattice class, shown in the listing below, inherits from Lattice and is used as a base class for lattices that are implemented in terms of one or more trees. Of course, it should be called TreeBasedLattice instead; the shorter name was probably chosen before we discovered the marvels of automatic completion—or English grammar.

    template <class Impl>
    class TreeLattice : public Lattice,
                        public CuriouslyRecurringTemplate<Impl> {
     public:
       TreeLattice(const TimeGrid& timeGrid, Size n);
       void initialize(DiscretizedAsset& asset, Time t) const {
           Size i = t_.index(t);
           asset.time() = t;
           asset.reset(this->impl().size(i));
       }
       void rollback(DiscretizedAsset& asset, Time to) const {
           partialRollback(asset,to);
           asset.adjustValues();
       }
       void partialRollback(DiscretizedAsset& asset, Time to) const {
           Integer iFrom = Integer(t_.index(asset.time()));
           Integer iTo = Integer(t_.index(to));
           for (Integer i=iFrom-1; i>=iTo; --i) {
               Array newValues(this->impl().size(i));
               this->impl().stepback(i, asset.values(), newValues);
               asset.time() = t_[i];
               asset.values() = newValues;
               if (i != iTo) // skip the very last adjustment
                   asset.adjustValues();
           }
       }
       void stepback(Size i, const Array& values,
                             Array& newValues) const {
           for (Size j=0; j<this->impl().size(i); j++) {
               Real value = 0.0;
               for (Size l=0; l<n_; l++) {
                   value += this->impl().probability(i,j,l) *
                            values[this->impl().descendant(i,j,l)];
               }
               value *= this->impl().discount(i,j);
               newValues[j] = value;
           }
       }
       Real presentValue(DiscretizedAsset& asset) const {
           Size i = t_.index(asset.time());
           return DotProduct(asset.values(), statePrices(i));
       }
       const Array& statePrices(Size i) const;
    };

This class template acts as an adapter between the Lattice class, from which it inherits the interface, and the Tree class template which will be used for the implementation. Once again, we used the Curiously Recurring Template Pattern (which wasn’t actually needed for trees, but it is in this case); the behavior of the lattice is written in terms of a number of methods that must be defined in derived classes. For greater generality, there is no mention of trees in the TreeLattice class. It’s up to derived classes to choose what kind of trees they should contain and how to use them.

The TreeLattice constructor is simple enough: it takes and stores the time grid and an integer n specifying the order of the tree (2 for binomial, 3 for trinomial and so on). It also performs a check or two, and initializes a couple of data members used for caching data; but I’ll gloss over that here.

The interesting part is the implementation of the Lattice interface, which follows the outline I gave back in the first post on this framework. The initialize method calculates the index of the passed time on the stored grid, sets the asset time, and finally passes the number of nodes on the corresponding tree level to the asset’s reset method. The number of nodes is obtained by calling the size method through CRTP; this is one of the methods that derived classes will have to implement, and (like all other such methods) has the same signature as the corresponding method in the tree classes.

The rollback and partialRollback methods perform the same work, with the only difference that rollback performs the adjustment at the final time and partialRollback doesn’t. Therefore, it’s only to be expected that the one is implemented in terms of the other; rollback performs a call to partialRollback, followed by another to the asset’s adjustValues method.

The rollback procedure is spelled out in partialRollback: it finds the indexes of the current time and of the target time on the grid, and it loops from one to the other. At each step, it calls the stepback method, which implements the actual numerical work of calculating the asset values on the i-th level of the tree from those on level i+1; then it updates the asset and, at all steps except the last, calls its adjustValues method.

The implementation of the stepback method defines, by using it, the interface that derived classes must implement (that’s currently the only sane way, since concepts were left out of C++11). It determines the value of the asset at each node by combining the values at each descendant node, weighed by the corresponding transition probability; the result is further adjusted by discounting it. All in all, the required interface includes the size method, which I’ve already shown; the probability and descendant methods, with the same signature as the tree methods of the same name; and the discount method, which takes the indexes of the desired level and node and returns the discount factor between that node and its descendants (assumed to be independent of the particular descendant).

Finally, the presentValue method is implemented by returning the dot-product of the asset values by the state prices at the current time on the grid. I’ll cheerfully ignore the way the state prices are calculated; suffice it to say that using them is somewhat more efficient than rolling the asset all the way back to \( t=0 \).

Now, why does the TreeLattice implementation calls methods with the same signature as those of a tree (thus forcing derived classes to define them) instead of just storing a tree and calling its methods directly? Well, that’s the straightforward thing to do if you have just an underlying tree; and in fact, most one-dimensional lattices will just forward the calls to the tree they store. However, it wouldn’t work for other lattices (say, two-dimensional ones); and in that case, the wrapper methods used in the implementation of stepback allow us to adapt the underlying structure, whatever that is, to their common interface.

The library contains instances of both kinds of lattices. Most—if not all—of those of the straightforward kind inherit from the TreeLattice1D class template, shown in the listing below.

    template <class Impl>
    class TreeLattice1D : public TreeLattice<Impl> {
      public:
        TreeLattice1D(const TimeGrid& timeGrid, Size n);
        Disposable<Array> grid(Time t) const;
    };

It doesn’t define any of the methods required by TreeLattice; and the method it does implement (the grid method, defined as pure virtual in the Lattice base class) actually requires another one, namely, the underlying method. All in all, this class does little besides providing a useful categorization; the storage and management of the underlying tree is, again, left to derived classes. (One might argue for including a default implementation of the required methods in TreeLattice1D. This would probably make sense; it would make it easier to implement derived classes in the most common cases, and could be overridden if a specific lattice needed it.)

One such class is the inner OneFactorModel::ShortRateTree class, shown in the next listing.

    class OneFactorModel::ShortRateTree
        : public TreeLattice1D<OneFactorModel::ShortRateTree> {
      public:
        ShortRateTree(
              const boost::shared_ptr<TrinomialTree>& tree,
              const boost::shared_ptr<ShortRateDynamics>& dynamics,
              const TimeGrid& timeGrid)
        : TreeLattice1D<OneFactorModel::ShortRateTree>(timeGrid, 3),
          tree_(tree), dynamics_(dynamics) {}
        Size size(Size i) const {
            return tree_->size(i);
        }
        Real underlying(Size i, Size index) const {
            return tree_->underlying(i, index);
        }
        Size descendant(Size i, Size index, Size branch) const {
            return tree_->descendant(i, index, branch);
        }
        Real probability(Size i, Size index, Size branch) const {
            return tree_->probability(i, index, branch);
        }
        DiscountFactor discount(Size i, Size index) const {
            Real x = tree_->underlying(i, index);
            Rate r = dynamics_->shortRate(timeGrid()[i], x);
            return std::exp(-r*timeGrid().dt(i));
        }
      private:
        boost::shared_ptr<TrinomialTree> tree_;
        boost::shared_ptr<ShortRateDynamics> dynamics_;
    };

Its constructor takes a trinomial tree, built by any specific short-rate model according to its dynamics; an instance of the ShortRateDynamics class, which I’ll gloss over; and a time grid, which could have been extracted from the tree so I can’t figure out why we pass it instead. The grid is passed to the base-class constructor, together with the order of the tree (which is 3, of course); the tree and the dynamics are stored as data members.

As is to be expected, most of the required interface is implemented by forwarding the call to the corresponding tree method. The only exception is the discount method, which doesn’t have a corresponding tree method; it is implemented by asking the tree for the value of its underlying value at the relevant node, by retrieving the short rate from the dynamics, and by calculating the corresponding discount factor between the time of the node and the next time on the grid.

Note that, by modifying the dynamics, it is possible to change the value of the short rate at each node while maintaining the structure of the tree unchanged. This is done in a few models in order to fit the tree to the current interest-rate term structure; the ShortRateTree class provides another constructor, not shown here, that takes additional parameters to perform the fitting procedure.

As an example of the second kind of lattice, have a look at the TreeLattice2D class template, shown in the listing below. It acts as base class for lattices with two underlying variables, and implements most of the methods required by TreeLattice. (In this, it differs from TreeLattice1D which didn’t implement any of them. We might have had a more symmetric hierarchy by leaving TreeLattice2D mostly empty and moving the implementation to a derived class. At this time, though, it would sound a bit like art for art’s sake.)

    template <class Impl, class T = TrinomialTree>
    class TreeLattice2D : public TreeLattice<Impl> {
      public:
        TreeLattice2D(const boost::shared_ptr<T>& tree1,
                      const boost::shared_ptr<T>& tree2,
                      Real correlation)
        : TreeLattice<Impl>(tree1->timeGrid(),
                            T::branches*T::branches),
          tree1_(tree1), tree2_(tree2), m_(T::branches,T::branches),
          rho_(std::fabs(correlation)) { ... }
        Size size(Size i) const {
            return tree1_->size(i)*tree2_->size(i);
        }
        Size descendant(Size i, Size index, Size branch) const {
            Size modulo = tree1_->size(i);

            Size index1 = index % modulo;
            Size index2 = index / modulo;
            Size branch1 = branch % T::branches;
            Size branch2 = branch / T::branches;

            modulo = tree1_->size(i+1);
            return tree1_->descendant(i, index1, branch1) +
                   tree2_->descendant(i, index2, branch2)*modulo;
        }
        Real probability(Size i, Size index, Size branch) const {
            Size modulo = tree1_->size(i);

            Size index1 = index % modulo;
            Size index2 = index / modulo;
            Size branch1 = branch % T::branches;
            Size branch2 = branch / T::branches;

            Real prob1 = tree1_->probability(i, index1, branch1);
            Real prob2 = tree2_->probability(i, index2, branch2);
            return prob1*prob2 + rho_*(m_[branch1][branch2])/36.0;
        }
      protected:
        boost::shared_ptr<T> tree1_, tree2_;
        Matrix m_;
        Real rho_;
    };

The two variables are modeled by correlating the respective trees. Now, I’m sure that any figure I might draw would only add to the confusion. However, the idea is that the state of the two variables is expressed by a pair of node from the respective trees; that the transitions to be considered are those from pair to pair; and that all the possibilities are enumerated so that they can be retrieved by means a single index and thus can match the required interface.

For instance, let’s take the case of two trinomial trees. Let’s say we’re at level \( i \) (the two trees must have the same time grid, or all bets are off). The first variable has a value that corresponds to node \( j \) on its tree, while the second sits on node \( k \). The structure of the first tree tells us that on next level, the first variable might go to nodes \( j’_0 \), \( j’_1 \) or \( j’_2 \) with different probabilities; the second tree gives us \( k’_0 \), \( k’_1 \) and \( k’_2 \) as target nodes for the second variable. Seen as transition between pairs, this means that we’re at \( (j,k) \) on the current level and that on the next level we might go to any of \( (j’_0,k’_0) \), \( (j’_1,k’_0) \), \( (j’_2,k’_0) \), \( (j’_0,k’_1) \), \( (j’_1,k’_1) \), and so on until \( (j’_2,k’_2) \) for a grand total of \( 3 \times 3 = 9 \) possibilities. By enumerating the pairs in lexicographic order like I just did, we can give \( (j’_0,k’_0) \) the index \( 0 \), \( (j’_1,k’_0) \) the index \( 1 \), and so on until we give the index \( 8 \) to \( (j’_2,k’_2) \). In the same way, if on a given level there are \( n \) nodes on the first tree and \( m \) on the second, we get \( n \times m \) pairs that, again, can be enumerated in lexicographic order: the pair \( (j,k) \) is given the index \( k \times n + j \).

At this point, the implementation starts making sense. The constructor of TreeLattice2D takes and stores the two underlying trees and the correlation between the two variables; the base-class TreeLattice constructor is passed the time grid, taken from the first tree, and the order of the lattice, which equals the product of the orders of the two trees; for two trinomial trees, this is \( 3 \times 3 = 9 \) as above. (The current implementation assumes that the two trees are of the same type, but it could easily be made to work with two trees of different orders.) The constructor also initializes a matrix m_ that will be used later on.

The size of the lattice at a given level is the product \( n \times m \) of the sizes of the two trees, which translates in a straightforward way into the implementation of the size method.

Things get more interesting with the two following methods. The descendant method takes the level i of the tree; the index of the lattice node, which is actually the index of a pair \( (j,k) \) among all those available; and the index of the branch to take, which by the same token is a pair of branches. The first thing it does is to extract the actual pairs. As I mentioned, the passed index equals \( k \times n + j \), which means that the two underlying indexes can be retrieved as index%n and index/n. The same holds for the two branches, with \( n \) being replaced by the order of the first tree. Having all the needed indexes and branches, the code calls the descendant method on the two trees, obtaining the indexes \( j’ \) and \( k’ \) of the descendant nodes; then it retrieves the size \( n’ \) of the first tree at the next level; and finally returns the combined index \( k’ \times n’ + j’ \).

Up to a certain point, the probability method performs the same calculations; that is, until it retrieves the two probabilities from the two underlying trees. If the two variables were not correlated, the probability for the transition would then be the product of the two probabilities. Since this is not the case, a correction term is added which depends from the passed correlation (of course) and also from the chosen branches. I’ll gloss on the value of the correction as I already did on several other formulas.

This completes the implementation. Even though it contains a lot more code than its one-dimensional counterpart, TreeLattice2D is still an incomplete class. Actual lattices will have to inherit from it, close the CRTP loop, and implement the missing discount method.

I’ll close this post by mentioning one feature that is currently missing from lattices, but would be nice to have. If you turn back to the implementation of the stepback method, you’ll notice that it assumes that the values we’re rolling back are cash values; that is, it always discounts them. It could also be useful to roll values back without discounting. For instance, the probability that an option be exercised could be calculated by rolling it back on the trees without discounting, while adjusting it to \( 1 \) on the nodes where the exercise condition holds.