Greetings.

This week, another post in a series on the new finite-difference framework. The end of the chapter doesn’t look so far away, now. And if you missed the series so far, the first post is here.

In other news, I’m just back from the QuantLib User Meeting in Düsseldorf. I’ll try to write a trip report shortly; in the meantime, you can check the QuantLib Twitter feed for pictures and a few links.

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 buttons for that are in the footer. Also, make sure to check my Training page.

Examples: Black-Scholes operators

Most full-featured operators in the library are not inherited directly from the FdmLinearOp class or from one of the basic operators. Instead, they inherit from the FdLinearOpComposite class (shown in the listing below), contain one or more basic operators as data members, and use them to implement their own behavior.

    class FdmLinearOpComposite : public FdmLinearOp {
      public:
        virtual Size size() const = 0;

        virtual void setTime(Time t1, Time t2) = 0;

        virtual Array apply_mixed(const Array& r) const = 0;
        virtual Array apply_direction(Size direction,
                                      const Array& r) const = 0;
        virtual Array solve_splitting(Size direction,
                                      const Array& r,
                                      Real s) const = 0;
        virtual Array preconditioner(const Array& r,
                                     Real s) const = 0;
    };

As you can see, the FdmLinearOpComposite class augments the interface of operators with a few more methods that derived classes must implement, the simplest being size (which must return the number of dimensions of the operator).

The other methods deserve a bit more attention. The setTime method implements the same idea I described in a previous post: when called, it should modify or rebuild a time-dependent operator so that its elements correspond to the correct time. However, there are a couple of enhancements. First, in the old framework the machinery was added to the TridiagonalOperator class, that is, to the basic building block for operators; in the new one, basic operators are constant (and simpler) and the setTime method was added to a higher-level interface. Second, the operator will be used over a given step, and this interface allows one to pass both its start and end time; this allows operators to provide a better discretization (for instance, by averaging or integrating) than simply the value at a single time.

The remaining methods are used to support Alternating Direction Implicit (ADI) schemes. I’m not going to describe them in any detail here (a summary is available, e.g., in [1], and you all know how to use Google anyway); the basic idea is that an operator \( A \) is decomposed as \( A = A_m + A_0 + A_1 + \ldots + A_{N-1} \), where \( A_m \) represents the mixed derivatives and each \( A_i \) represents the derivatives along the \( i \)-th direction, and those components are used separately. Instead of having the operator create the actual components and let the schemes apply them, the interface declares corresponding methods: thus, the apply_mixed method will apply the \( A_m \) component and the apply_direction method will apply one of the \( A_i \) depending on the direction passed to the method. The solve_splitting and preconditioner methods are also used in ADI schemes.

As a first example, you can look at the FdmBlackScholesOp class, shown in the next listing.

    class FdmBlackScholesOp : public FdmLinearOpComposite {
      public:
        FdmBlackScholesOp(
            const shared_ptr<FdmMesher>& mesher,
            const shared_ptr<GeneralizedBlackScholesProcess>&,
            Real strike,
            bool localVol = false,
            Real illegalLocalVolOverwrite = -Null<Real>(),
            Size direction = 0);

        Size size() const { return 1; }
        void setTime(Time t1, Time t2) {
            if (localVol_) {
                /* `more calculations` */
            } else {
                Rate r = rTS_->forwardRate(t1, t2, Continuous);
                Rate q = qTS_->forwardRate(t1, t2, Continuous);
                Real v = volTS_->blackForwardVariance(
                                     t1, t2, strike_)/(t2-t1);
                mapT_ = /* `put them together` */;
            }
        }

        Array apply(const Array& r) const {
            return mapT_.apply(r);
        }
        Array apply_mixed(const Array& r) const {
            return Array(r.size(), 0.0);
        }
        Array apply_direction(Size direction,
                              const Array& r) const {
            return (direction == direction_) ? mapT_.apply(r)
                                             : Array(r.size(), 0.0);
        }
        Array solve_splitting(Size direction,
                              const Array& r, Real s) const;
        Array preconditioner(const Array& r, Real s) const;
      private:
        TripleBandLinearOp mapT_;
        /* ... other data members ... */
    };

As expected, it inherits from FdmLinearOpComposite and implements its required interface. The constructor takes quite a few parameters and stores them in the corresponding data members for later use; among them, it takes a given direction, which let us specify the axis along which this operator works and thus allows it to be used as a building block in other operators.

The underlying differential operator, stored in the mapT_ data member, is built inside the setTime method. As I mentioned, having both the start time t1 and the end time t2 of the step allows the code to provide a more accurate discretization; namely, it can ask the term structures for the exact forward rates and variance between those two times instead of picking an instantaneous value. Depending on whether we want to use local volatility, the calculation of the diffusion term differs; but in both cases we end up building the Black-Scholes operator

\[L_{BS}(t) = - \left(r - q - \frac{\sigma^2}{2} \right) \frac{\partial}{\partial x} - \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2} + r\]

and storing it into mapT_ so that it can be used elsewhere.

Once the final operator is built, the implementation of the remaining methods is straightforward enough. The apply method applies mapT_ to the passed array and return the results. Since this is a one-dimensional operator, there are no cross-terms, therefore apply_mixed returns a null array. Finally, the apply_direction method applies mapT_ to the passed array when the direction equals the one specified in the constructor (because that’s the one direction along which the entire operator works) and instead returns a null array when the direction is different (because the operator has no corresponding component). The implementations of the other methods work in a similar way.

As a second example, the Fdm2dBlackScholesOp class, shown in the listing below, builds a two-dimensional Black-Scholes operator by combining a pair of one-dimensional operators with their correlation.

    class Fdm2dBlackScholesOp : public FdmLinearOpComposite {
      public:
        Fdm2dBlackScholesOp(
            const shared_ptr<FdmMesher>& mesher,
            const shared_ptr<GeneralizedBlackScholesProcess>& p1,
            const shared_ptr<GeneralizedBlackScholesProcess>& p2,
            Real correlation,
            Time maturity,
            bool localVol = false,
            Real illegalLocalVolOverwrite = -Null<Real>());
           
        Size size() const { return 2; }
        void setTime(Time t1, Time t2) {
            opX_.setTime(t1, t2);
            opY_.setTime(t1, t2);

            corrMapT_ = /* ... other calculations ... */
        }

        Array apply(const Array& x) const {
            return opX_.apply(x) + opY_.apply(x) + apply_mixed(x);
        }
        Array apply_mixed(const Array& x) const {
            return corrMapT_.apply(x) + currentForwardRate_*x;
        }
        Array apply_direction(Size direction,const Array& x) const {
            if (direction == 0)
                return opX_.apply(x);
            else if (direction == 1)
                return opY_.apply(x);
            else
                QL_FAIL("direction is too large");
        }
        Array solve_splitting(Size direction,
                              const Array& x, Real s) const;
        Array preconditioner(const Array& r, Real s) const;
      private:
        FdmBlackScholesOp opX_, opY_;
        NinePointLinearOp corrMapT_;
        /* ... other data members ... */
    };

Again, the constructor stores the passed data, while the setTime method builds the operators; directly in the case of the correlation corrMapT_, and forwarding to the corresponding methods in the case of the two one-dimensional operators opX_ and opY_.

The other methods are, again, simple. By linearity of the operators, the apply method works by applying each component to the passed array in turn and returning the sum of the results; apply_mixed uses the mixed operator corrMapT_ plus a constant term; and apply_direction selects and applies the correct one-dimensional component.

If you’re interested, the library provides other operators you can examine. Most of them turn out to have the same structure: for instance, the FdmHestonOp class, which implements a helper operator for each underlying variable and puts them together in the final operator, or the FdmHestonHullWhiteOp, that you can inspect as an example of three-dimensional operator.

Bibliography

[1] C. S. L. de Graaf. Finite Difference Methods in Derivatives Pricing under Stochastic Volatility Models. Master’s thesis, Mathematisch Instituut, Universiteit Leiden, 2012.