Welcome back.

Spring is well underway, and my young colleagues’ fancy lightly turns to the dating apps on their phones. I’m married with children, so I finished my series on the new finite-difference framework instead. The collected series forms chapter 8 of Implementing QuantLib, whose main content is now finished. I’m still missing about half of the appendix, after which I’ll have the first complete version. But of course, Leanpub being Leanpub, you can buy the book now and get further versions for free when I update it. Neat, huh?

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.

Schemes and solvers

As in the old framework, schemes provide a discretization of the time derivative in the equation we want to solve; you can go back to an earlier post for more details. The new framework provides reimplementations of schemes used in the old one, such as explicit and implicit Euler, as well as a number of new ones I won’t describe; as usual, you can look them up in the library.

Like in the old framework, they have no common base class; they all courteously agree to implement step and setStep methods with the same semantics described in the section I just quoted, but there’s no inheritance relationship holding them to any contract. This is a bit of a problem, given that the new framework usually eschews generic programming and templates in favor of more object-oriented techniques; the absence of a base class, of course, prevents us from writing non-template code taking an instance of a generic scheme.

The way the framework works around this is to define the FdmSchemeDesc class, shown in the following listing. Its instances can be passed around instead of the scheme they describe, and contain the type of the chosen scheme and the parameters that might be needed to instantiate it. The available types are enumerated in the class declaration, and yes, that’s not optimal. I’ll get back to this later.

    struct FdmSchemeDesc {
        enum FdmSchemeType { HundsdorferType, DouglasType,
                             CraigSneydType, ModifiedCraigSneydType,
                             ImplicitEulerType, ExplicitEulerType };

        FdmSchemeDesc(FdmSchemeType type, Real theta, Real mu);

        const FdmSchemeType type;
        const Real theta, mu;

        static FdmSchemeDesc Douglas();
        static FdmSchemeDesc ImplicitEuler();
        static FdmSchemeDesc ExplicitEuler();
        static FdmSchemeDesc CraigSneyd();
        static FdmSchemeDesc ModifiedCraigSneyd();
        static FdmSchemeDesc Hundsdorfer();
        static FdmSchemeDesc ModifiedHundsdorfer();
    };

Besides the enumeration, the class declares a bunch of static methods returning a few pre-built schemes with sensible parameters; if you want a custom one instead, you can build an instance with any old values.

Note also that this class and the others described in this section are not templates; they’re written to work specifically with the FdmLinearOpComposite class. My guess is that their authors learned the lesson from the old framework, which was written for a generic operator and turned out to be only ever used with the TridiagonalOperator class.

Now, if the type enumeration seemed a bit icky to you, you might want to brace yourself for the next listing and the FdmBackwardSolver class.

    class FdmBackwardSolver {
      public:
        FdmBackwardSolver(
            const shared_ptr<FdmLinearOpComposite>& map,
            const FdmBoundaryConditionSet& bcSet,
            const shared_ptr<FdmStepConditionComposite> condition,
            const FdmSchemeDesc& schemeDesc);

        void rollback(array_type& rhs,
                      Time from, Time to,
                      Size steps, Size dampingSteps) {
            const Time deltaT = from - to;
            const Size allSteps = steps + dampingSteps;
            const Time dampingTo =
                from - (deltaT*dampingSteps)/allSteps;

            if (dampingSteps && schemeDesc_.type !=
                                FdmSchemeDesc::ImplicitEulerType) {
                ImplicitEulerScheme implicitEvolver(map_, bcSet_);
                FiniteDifferenceModel<ImplicitEulerScheme>
                dampingModel(implicitEvolver,
                             condition_->stoppingTimes());
                dampingModel.rollback(rhs, from, dampingTo,
                                      dampingSteps, *condition_);
            }

            switch (schemeDesc_.type) {
              case FdmSchemeDesc::HundsdorferType:
                {
                    HundsdorferScheme hsEvolver(schemeDesc_.theta,
                                                schemeDesc_.mu,
                                                map_, bcSet_);
                    FiniteDifferenceModel<HundsdorferScheme>
                    hsModel(hsEvolver, condition_->stoppingTimes());
                    hsModel.rollback(rhs, dampingTo, to,
                                     steps, *condition_);
                }
                break;
              // other cases, not shown
            }
        }
    };

Its constructor collects and stores an operator, a set of boundary and step conditions, and a scheme description. They are used in the rollback method, which puts them together in a finite-difference model, takes an array of initial values, and runs the model between the two given times from and to in the given number of steps, possibly with a few initial damping steps. And, as you can see, the implementation of rollback is a big switch on type.

Was that a gasp? Yes, I know. Both the FdmSchemeDesc class and the rollback method need to be changed if we want to add a new scheme, which is a violation of the open-closed principle. It’s not like the authors of the code were sloppy, though. Let’s say we go ahead and inherit all schemes from a base class, so we can use a polymorphic pattern such as Strategy to select one. There are two forces at play here. First, if the user must be able to choose it, the scheme instance would have to be passed to the engine from outside. Second, the operator and the boundary conditions, which the scheme requires, are created inside the engine, because having the user write code to create them would cause duplication, and thus the scheme must be created inside the engine, too. Tricky, isn’t it?

Two possible solutions come to mind. One is to use some kind of factory: the scheme might have, say, a clone method that returns another instance with the same type and parameters and with the passed boundary and step conditions added (we could also have factory classes, but that would mean a parallel hierarchy, which is never good). Another is to change the interface of the classes we’ve seen, so that the scheme doesn’t contain the conditions; since backward compatibility prevents us from changing the FiniteDifferenceModel class, this requires writing a new model class, too. Both would have been possible, and the maintainer who reviewed and merged the framework when it was contributed should have seen the problem and should have suggested one of them. That maintainer was me, in case you’re curious. Live and learn, I guess.

And after the FdmBackwardSolver class, we’re on the home stretch; the remaining classes add a couple of layers that make it more convenient to use the solver, but not a lot of new logic. To the more observant among you: yes, in this section I abandoned my initial intent to write a top-down description of the components. I think it makes more sense this way.

First, the FdmSolverDesc structure shown in the listing that follows groups together the arguments needed to initialize a solver so that we can pass them around together. The operator is not included; that’s because, as you’ll see in a minute, we’ll define higher-level solvers that define the operator and can be reused for different instruments by passing them different FdmSolverDesc instances.

    struct FdmSolverDesc {
        const boost::shared_ptr<FdmMesher> mesher;
        const FdmBoundaryConditionSet bcSet;
        const boost::shared_ptr<FdmStepConditionComposite> condition;
        const boost::shared_ptr<FdmInnerValueCalculator> calculator;
        const Time maturity;
        const Size timeSteps;
        const Size dampingSteps;
    };

Second, the Fdm1DimSolver class sketched in the next listing takes care of some of the plumbing around a solver.

    class Fdm1DimSolver : public LazyObject {
      public:
        Fdm1DimSolver(const FdmSolverDesc& solverDesc,
                      const FdmSchemeDesc& schemeDesc,
                      const shared_ptr<FdmLinearOpComposite>& op)
        : /* ... */ {
            // ... get mesher, calculator etc. from solver description ...
            FdmLinearOpIterator end = layout->end();
            for (FdmLinearOpIterator i = layout->begin();
                                     i != end; ++i) {
                initialValues_[i.index()] =
                    calculator->avgInnerValue(i, maturity);
                x_[i.index()] = mesher->location(i, 0);
            }
        }

        Real interpolateAt(Real x) const {
            calculate();
            return (*interpolation_)(x);
        }
        // ...other inspectors for derivatives
      protected:
        void performCalculations() const {
            Array rhs(initialValues_.size());
            std::copy(initialValues_.begin(), initialValues_.end(),
                      rhs.begin());

            FdmBackwardSolver(op_, solverDesc_.bcSet,
                              conditions_, schemeDesc_)
                .rollback(rhs, solverDesc_.maturity, 0.0,
                          solverDesc_.timeSteps,
                          solverDesc_.dampingSteps);

            std::copy(rhs.begin(), rhs.end(), resultValues_.begin());
            interpolation_ = shared_ptr<CubicInterpolation>(
                new MonotonicCubicNaturalSpline(
                    x_.begin(), x_.end(), resultValues_.begin()));
        }
    };

Its constructor takes a solver description, a scheme description and an operator, stores them, and prepares a number of other data members for the calculations. In particular, it uses the layout in the solver description to allocate the correct size for the initial values, and the calculator to fill them; also, it extracts from the mesher the values of the underlying variable on the chosen grid.

The calculations are implemented inside the performCalculations method (yes, the class inherits from LazyObject, which puzzles me; I’ll come back to this). As I mentioned, it’s mostly plumbing: the code makes a copy of the initial values, creates an instance of FdmBackwardSolver, uses it to roll back the values, stores the results, and wraps them in an interpolation. Methods such as interpolateAt make the results available to client code; other inspectors, not shown, give access to the first and second derivative of the interpolation, as well as the time derivative of the result.

Going back one step: I gave the matter some thought, but I don’t see the need for inheriting from LazyObject here. For one thing, the constructor of Fdm1DimSolver doesn’t register with anything, which kind of defeats the purpose of the pattern (registration might happen in derived classes, but a quick search in the library shows that Fdm1DimSolver doesn’t have any). Moreover, we’ll see that solvers are normally created by an engine during calculation, used right there, and discarded afterwards. All in all, we could have put the whole calculation in the constructor, or even turned the class into a function returning some kind of result structure containing the interpolation and some inspectors to access it.

I should mention that Fdm1DimSolver has siblings: the framework also defines the Fdm2DimSolver and Fdm3DimSolver classes, which perform the same calculations and produce a 2-D or 3-D interpolation, respectively. The interface is almost the same, with some difference due to the higher dimensions; for instance, the interpolateAt method takes as arguments x and y in the 2-D case and x, y and z in the 3-D case.

Finally, specific solvers are tasked with creating an operator and using it with one of the generic solvers I described. It’s the case of the FdmBlackScholesSolver class, shown in the next listing.

    class FdmBlackScholesSolver : public LazyObject {
      public:
        FdmBlackScholesSolver(
            const Handle<GeneralizedBlackScholesProcess>& process,
            Real strike,
            const FdmSolverDesc& solverDesc,
            const FdmSchemeDesc& schemeDesc,
            bool localVol = false,
            Real illegalLocalVolOverwrite = -Null<Real>())
        : /* ... */ {
            registerWith(process_);
        }

        Real valueAt(Real s) const {
            calculate();
            return solver_->interpolateAt(std::log(s));
        }
        Real deltaAt(Real s) const {
            calculate();
            return solver_->derivativeX(std::log(s))/s;
        }
        Real gammaAt(Real s) const;
        Real thetaAt(Real s) const;

      protected:
        void performCalculations() const {
            const boost::shared_ptr<FdmBlackScholesOp> op(
                new FdmBlackScholesOp(
                    solverDesc_.mesher, process_.currentLink(),
                    strike_, localVol_, illegalLocalVolOverwrite_));

            solver_ = boost::shared_ptr<Fdm1DimSolver>(
                new Fdm1DimSolver(solverDesc_, schemeDesc_, op));
        }

      private:
        // other data members, not shown
        mutable boost::shared_ptr<Fdm1DimSolver> solver_;
    };

Besides the usual descriptions and a few other parameters for the model, it takes a handle to a Black-Scholes process that it stores and with which it registers for notifications; again, this would make more sense if engines kept the solver around instead of destroying it. Its performCalculations method creates an operator based on the process and uses it to build a 1-D solver, which is then stored. Inspectors such as valueAt and deltaAt call calculate first, which triggers the construction of the 1-D solver, and then delegate to the corresponding methods of the stored solver, which in turn trigger the rollback calculation inside it. The inspectors also take care of the transformation between the underlying value, used in the interface, and its logarithm, used in the process.

As I mentioned earlier, the solver takes the solver description from the outside and thus can be used for different engines, provided that the underlying of the option follows the same dynamics. In the listing in this post, you can see it used for a vanilla put or call option; by changing the boundary conditions, it can be used, and in fact it is used, for barrier options; and by changing the initial conditions, it might be used for digital options or other kinds of instruments. You can see how this design is more flexible than the old, inheritance-based one; or at least, that it is much simpler to get flexibility.

Multi-dimensional solvers, which I won’t describe, work in the same way. If you’re interested, you can look at the FdmBatesSolver class for a 2-D example and at the FdmHestonHullWhiteSolver class for a 3-D one. And, as the saying goes, there’s plenty more where those come from.

With solvers, we finally have all the pieces we need to build our engines—and round up this chapter, I might add. Looking back at the listing of the engine we took as example, we can now see how a finite-difference engine works. It takes the scheme and the other required parameters at construction, and then performs its calculations by building the required mesher, the initial, step and boundary conditions, and finally the solver from which the results can be extracted.