Welcome back.

This post is the third in a series (started here) that covers chapter 6 of my book. This time, an example: the Black-Scholes process (in which I’ll point out all the flaws in our code, thus showing that I’m not cut out for marketing).

You still have a few days (until the end of this month) to get an early-bird discount for my Introduction to QuantLib Development course, which will be on September 22nd to 24th in London. I’ve talked about it at length in a previous post, so I’ll shut up now.

Follow me on Twitter if you want to be notified of new posts, or add me to your circles, or subscribe via RSS: the buttons for that are in the footer. Also, make sure to check my Training page.

Example: the Black-Scholes process

And now, an example or three. Not surprisingly, the first is the one-dimensional Black-Scholes process. You’ve already met briefly the corresponding class (by the unwieldy name of GeneralizedBlackScholesProcess) in this post, where it was passed as an argument to a pricing engine for a European option; I’ll go back to that code later on, when I discuss a few shortcomings of the current design and possible future solutions. For the time being, let’s have a look at the current implementation, sketched in listing 6.6.

Listing 6.6: Partial implementation of the GeneralizedBlackScholesProcess class.

    class GeneralizedBlackScholesProcess
        : public StochasticProcess1D {
            const Handle<Quote>& x0,
            const Handle<YieldTermStructure>& dividendTS,
            const Handle<YieldTermStructure>& riskFreeTS,
            const Handle<BlackVolTermStructure>& blackVolTS,
            const boost::shared_ptr<discretization>& d =
        Real x0() const {
            return x0_->value();
        Real drift(Time t, Real x) const {
            Real sigma = diffusion(t,x);
            Time t1 = t + 0.0001;
            return riskFreeRate_->forwardRate(t,t1,Continuous,...)
                 - dividendYield_->forwardRate(t,t1,Continuous,...)
                 - 0.5 * sigma * sigma;
        Real diffusion(Time t, Real x) const;
        Real apply(Real x0, Real dx) const {
            return x0 * std::exp(dx);
        Real expectation(Time t0, Real x0, Time dt) const {
            QL_FAIL("not implemented");
        Real evolve(Time t0, Real x0, Time dt, Real dw) const {
            return apply(x0, discretization_->drift(*this,t0,x0,dt) +
        const Handle<Quote>& stateVariable() const;
        const Handle<YieldTermStructure>& dividendYield() const;
        const Handle<YieldTermStructure>& riskFreeRate() const;
        const Handle<BlackVolTermStructure>& blackVolatility() const;
        const Handle<LocalVolTermStructure>& localVolatility() const;

Well, first of all, I should probably explain the length of the class name. The “generalized” bit doesn’t refer to some extension of the model, but to the fact that this class was meant to cover different specific processes (such as the Black-Scholes process proper, the Black-Scholes-Merton process, the Black process and so on); the shorter names were reserved for the specializations.

The constructor of this class takes handles to the full set of arguments needed for the generic formulation: a term structure for the risk-free rate, a quote for the varying market value of the underlying, another term structure for its dividend yield, and a term structure for its implied Black volatility. In order to implement the StochasticProcess interface, the constructor also takes a discretization instance. All finance-related inputs are stored in the constructed class instance and can be retrieved by means of inspectors.

You might have noticed that the constructor takes a Black volatility \( \sigma_B(t,k) \), which is not what is needed to return the \( \sigma(t,x) \) diffusion term. The process performs the required conversion internally, using different algorithms depending on the type of the passed Black volatility (constant, depending on time only, or with a smile). The type is checked at run-time with a dynamic cast, which is not pretty but in this case is simpler than a Visitor pattern. At this time, the conversion algorithms are built into the process and can’t be changed by the user; nor it is possible to provide a precomputed local volatility. The class uses the Observer pattern to determine when to convert the Black volatility into the local one. Recalculation is lazy; meaning that it happens only when the local volatility is actually used and not immediately upon a notification from an observable.

The next methods, drift and diffusion, show the crux of this class: it’s actually a process for the logarithm of the underlying that masquerades as a process for the underlying itself. Accordingly, the term returned by drift is the one for the logarithm, namely, \( r(t) - q(t) - \frac{1}{2}\sigma^2(t,x) \), where the rates and the local volatility are retrieved from the corresponding term structures (this is another problem in itself, to which I’ll return). The diffusion method returns the local volatility.

The choice of \( \log(x) \) as the actual stochastic variable brings quite a few consequences. First of all, the variable we’re using is not the one we’re exposing. The quote we pass to the constructor holds the market value of the underlying; the same value is returned by the x0 method; and the values passed to and returned from other methods such as evolve are for the same quantity. In retrospect, this was a bad idea (yes, I can hear you say “D’oh.”) We were led to it by the design we had developed, but it should have been a hint that some piece was missing.

Other consequences can be seen all over the listing. The choice of \( \log(x) \) as the working variable is the reason why the apply method is virtual, and indeed why it exists at all: applying to \( x \) a drift \( \Delta \) meant for \( log(x) \) is not done by adding them, but by returning \( x\exp(\Delta) \). The expectation method is affected, too: its default implementation would apply the drift as above, thus returning \( \exp(E[\log(x)]) \) instead of \( E[x] \). To prevent it from returning the wrong value, the method is currently disabled (it raises an exception). In turn, the evolve method, whose implementation relied on expectation, had to be overridden to work around it. (The evolve method might have been overridden for efficiency, even if expectation were available. The default implementation in StochasticProcess1D would return \( x\exp(\mu dt)\exp(\sigma dw) \), while the overridden one calculates \( x\exp(\mu dt + \sigma dw) \).)

On top of the above smells, we have a performance problem—which is known to all those that tried a Monte Carlo engine from the library, as well as to all the people on the Wilmott forums to whom the results were colorfully described. Using apply, the evolve method performs an exponentiation at each step; but above all, the drift and diffusion methods repeatedly ask term structures for values. If you remember this post, this means going through at least a couple of levels of virtual method calls, sometimes (if we’re particularly unlucky) retrieving a rate from discount factors that were obtained from the same rate to begin with.

Finally, there’s another design problem. As I mentioned previously, we already used the Black-Scholes process in the AnalyticEuropeanEngine shown as an example in this post. However, if you look at the engine code in the library, you’ll see that it doesn’t use the process by means of the StochasticProcess interface; it just uses it as a bag of quotes and term structures. This suggests that the process is a jack of too many trades.

How can we fix it, then? We’ll probably need a redesign. It’s all hypothetical, of course; there’s no actual code yet, just a few thoughts I had while I was reviewing the current code for this chapter. But it could go as follows.

First of all (and at the risk of being told once again that I’m overengineering) I’d separate the stochastic-process part of the class from the bag-of-term-structures part. For the second one, I’d create a new class, probably called something like BlackScholesModel. Instances of this class would be the ones passed to pricing engines, and would contain the full set of quotes and term structures; while we’re at it, the interface should also allow the user to provide a local volatility, or to specify new ways to convert the Black volatility with some kind of Strategy pattern.

From the model, we could build processes. Depending on the level of coupling and the number of classes we prefer, we might write factory classes that take a model and return processes; we might add factory methods to the model class that return processes; or we might have the process constructors take a model. In any case, this could allow us to write process implementations that could be optimized for the given simulation. For instance, the factory (whatever it might be) might take as input the time nodes of the simulation and return a process that precomputed the rates \( r(t_i) \) and \( q(t_i) \) at those times, since they don’t depend on the underlying value; if a type check told the factory that the volatility doesn’t have smile, the process could precompute the \( \sigma(t_i) \), too.

As for the \( x \) vs \( \log(x) \) problem, I’d probably rewrite the process so that it’s explicitly a process for the logarithm of the underlying: the x0 method would return a logarithm, and so would the other methods such as expectation or evolve. The process would gain in consistency and clarity; however, since it would be no longer guaranteed that a process generates paths for the underlying, client code taking a generic process would also need a function or a function object that converts the value of the process variable into a value of the underlying. Such function might be added to the StochasticProcess interface, or passed to client code along with the process.