Hello everybody.

This week, I finally continue the series of posts on the new finite-difference framework. It seems we’re on the home stretch.

After you read it, or before, or even in the meantime, you might download QuantLib 1.9.2. It’s exactly the same as QuantLib 1.9.1, except that the tests aren’t broken. For some reason, that seemed to confuse people.

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.

Initial conditions

In the old framework, as you might remember, base engines declared a virtual method that derived classes had to implement and which applied the initial condition of the problem. In the new framework, in keeping with its plug-and-play nature, calculating the initial condition is left instead to separate objects that can be reused in different engines. The base class for such objects is FdmInnerValueCalculator, shown in the listing below.

    class FdmInnerValueCalculator {
        virtual ~FdmInnerValueCalculator() {}

        virtual Real innerValue(const FdmLinearOpIterator& iter,
                                Time t) = 0;
        virtual Real avgInnerValue(const FdmLinearOpIterator& iter,
                                   Time t) = 0;

Its interface contains the innerValue method, which returns the initial value at a given point in the mesh (denoted by an iterator) and the avgInnerValue method, which returns an average of the initial value around a given point and thus might provide a way to smooth discontinuities in the payoff. They’re both declared as pure virtual; this makes complete sense for innerValue, and is not wrong for avgInnerValue either. In hindsight, though, a few derived classes give up the calculation of the average as too difficult and return the value at the central point, i.e., the return value of innerValue, as an approximation; this might have been a default implementation. (But you could also argue that the approximation should be a deliberate choice of whoever writes the derived class, and you wouldn’t be wrong either.)

The next listing contains an examples of concrete class implementing the interface. The FdmLogInnerValue class is used to calculate a payoff, depending on a single underlying variable, on a logarithmic mesh.

    class FdmLogInnerValue : public FdmInnerValueCalculator {
        FdmLogInnerValue(const shared_ptr<Payoff>& payoff,
                         const shared_ptr<FdmMesher>& mesher,
                         Size direction);

        Real innerValue(const FdmLinearOpIterator& iter, Time) {
            Real s = std::exp(mesher_->location(iter, direction_));
            return (*payoff_)(s);
        Real avgInnerValue(const FdmLinearOpIterator& iter, Time);

It takes the payoff itself, the mesh, and the index of the axis along which the underlying varies; the latter is needed because the mesh might be multi-dimensional. For instance, the underlying might be a stock value and the problem might model the stock and its volatility, or the stock and the interest rate, or all three. There’s no guarantee that the stock varies along the first axis.

As you probably expect, the innerValue method extracts from the mesh the value of the underlying variable corresponding to the passed iterator, and then passes it to the payoff. The avgInnerValue method does something similar, but with a combination of numerical integrals and boost::bind that I won’t bother to show here (and that I hope to eradicate when we can use lambdas instead).

A second example, FdmLogBasketInnerValue, is used for a payoff depending on the values of multiple underlying variables (and not on the total basket value, which probably makes the class name a misnomer).

    class FdmLogBasketInnerValue : public FdmInnerValueCalculator {
                const shared_ptr<BasketPayoff>& payoff,
                const shared_ptr<FdmMesher>& mesher);

        Real innerValue(const FdmLinearOpIterator& iter, Time) {
            Array x(mesher_->layout()->dim().size());
            for (Size i=0; i < x.size(); ++i)
                x[i] = std::exp(mesher_->location(iter, i));
            return (*payoff_)(x);
        Real avgInnerValue(const FdmLinearOpIterator& iter, Time) {
            return innerValue(iter, t);

Its innerValue method collects the underlying values from the several axis of the mesh—which, in this case, is assumed not to model any other variable—and passes them to the payoff. The avgInnerValue method bails out (poor thing) and calls innerValue.

Aside: dispelling magic.

If you look back at the creation of the inner-value calculator in the engine code shown in an earlier post, you’ll see that we pass the index of the dimension as the magic number 0. This, of course, is not ideal. The problem is that the info is hard-coded in the definition of the mesher, but it’s not accessible by code, and it’s not easy to think of a meaningful common interface that meshers might provide to retrieve it.

I don’t know if there’s a solution, other than the time-honored mitigation method of giving magic numbers a name; for instance, the FdmBlackScholesMesher class could provide a specific static variable such as x_dimension. However, good names for these variables might be hard to define; and the idea runs counter to a possible simplification I mentioned in yet another post, in which meshers could be non-polymorphic and derived classes would be replaced by factory functions.

Boundary conditions

The new framework doesn’t declare any new interface for boundary conditions; the FdmBoundaryConditionSet seen in the engine is defined as

    typedef OperatorTraits<FdmLinearOp>::bc_set

where OperatorTraits comes from the old framework and thus makes use of the BoundaryCondition class template I described back in the day; the only thing new is that we’re instantiating it with the FdmLinearOp class.

Of course, specific boundary conditions are implemented in this framework so that they work with the new classes. For instance, the next listing shows the interface of the FdmDirichletBoundary class. Its constructor takes the mesher that defines the domain of the problem, the value of the instrument on the boundary, the index of the axis whose boundary we want to control and an enumeration to specify the lower or upper side of the mesh. The implementation, not shown, also uses the mesher to identify the array elements to condition.

    class FdmDirichletBoundary
        : public BoundaryCondition<FdmLinearOp> {
        FdmDirichletBoundary(const shared_ptr<FdmMesher>& mesher,
                             Real valueOnBoundary,
                             Size direction, Side side);
        void applyBeforeApplying(operator_type&) const;
        void applyBeforeSolving(operator_type&, array_type&) const;
        void applyAfterApplying(array_type&) const;
        void applyAfterSolving(array_type&) const;
        void setTime(Time) {}

Step conditions

Finally, step conditions also inherit from a class template defined in the old framework (StepCondition, described here); and of course, they also use the classes defined in the new framework. One such condition is implemented in the FdmAmericanStepCondition class, shown in the listing below, which determines if an American option should be exercised.

    class FdmAmericanStepCondition : public StepCondition<Array> {
            const shared_ptr<FdmMesher>&,
            const shared_ptr<FdmInnerValueCalculator>&);

        void applyTo(Array& a, Time) const {
            shared_ptr<FdmLinearOpLayout> layout = mesher_->layout();

            for (FdmLinearOpIterator iter = layout->begin();
                                     iter != layout->end(); ++iter) {
                Real innerValue = calculator_->innerValue(iter, t);
                if (innerValue > a[iter.index()]) {
                    a[iter.index()] = innerValue;
        shared_ptr<FdmMesher> mesher_;
        shared_ptr<FdmInnerValueCalculator> calculator_;

The implementation is how you would expect: the applyTo method iterates over the mesh, and at each point it compares the intrinsic value given by the option payoff with the current value of the option and it updates the corresponding array element. The constructor takes and stores the information that applyTo needs.

The framework also defines the FdmStepConditionComposite class. It is a convenience class, which uses the Composite pattern to group multiple step conditions and pass them around as a single entity. Its interface is shown in the listing that follows; it defines the required applyTo method, which applies all conditions in sequence, as well as a number of utility methods. Most of them are generic enough; there are inspectors for the stored information, as well as a method to join multiple conditions (or groups thereof).

    class FdmStepConditionComposite : public StepCondition<Array> {
        typedef list<shared_ptr<StepCondition<Array> > > Conditions;

            const list<vector<Time> > & stoppingTimes,
            const Conditions & conditions);

        void applyTo(Array& a, Time t) const;

        const vector<Time>& stoppingTimes() const;
        const Conditions& conditions() const;

        static shared_ptr<FdmStepConditionComposite> joinConditions(
                    const shared_ptr<FdmSnapshotCondition>& c1,
                    const shared_ptr<FdmStepConditionComposite>& c2);

        static shared_ptr<FdmStepConditionComposite>
             const DividendSchedule& schedule,
             const shared_ptr<Exercise>& exercise,
             const shared_ptr<FdmMesher>& mesher,
             const shared_ptr<FdmInnerValueCalculator>& calculator,
             const Date& refDate,
             const DayCounter& dayCounter);

However, the vanillaComposite static method is a bit different. It is a factory method that takes a bunch of arguments and makes a ready-to-go composite condition that can be used in different problems; which is good, but it’s more specific than the rest of the class. Not that specific, mind you, since a quick search shows that it’s used in engines for a few instruments as different as barrier options and swaptions; but still, if we were to write more of them, I wouldn’t like them to proliferate inside this class. It might be cleaner to define them as free functions.

And we’re done for this week. In next post in the series, we’ll look at finite-difference schemes.