## Thursday, December 17, 2015

### Christmas break

Welcome back.

Just a short note to tell that, as the song says, I'll be home for Christmas. And unlike the song, it won't be only in my dreams.

Blogging will resume after the holidays (not that I have a very regular schedule anyway, but still).

Have a wonderful time, everybody.

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

## Thursday, December 10, 2015

### Screencast: my talk at the QuantLib user meeting 2015

Welcome back.

In last post, I promised more info on my talk on QuantLib, IPython Notebook and Docker at the QuantLib user meeting. Well, I managed to record a screencast on my laptop while I was presenting, so you can now watch and listen to it. It's just the good parts, too, as you can see the screen but not my face.

And by the way, most speakers have now contributed their slides; so when you're done here, go to the meeting page and have a good look.

Enjoy.

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

## Thursday, December 3, 2015

### Report from the QuantLib user meeting in Düsseldorf

Hello everybody.

I'm back from the QuantLib user meeting (it was last Monday and Tuesday). While I was there, I was told that it must be nice to see all the things that people did with the library.

Oh, yes. Definitely.

As usual, it was great to be there. A good variety of talks, good organization, and on top of that, the chance to see people that I usually just exchange emails with. As usual, a big thanks goes to Michael von der Driesch for preparing the meeting and for keeping it running smoothly, as well as this year's sponsors IKB and CompatibL.

Here is a short summary of the event. For more details, you can check the meeting page on the QuantLib site; slides are being collected and published, so you should find them all there in a week or two.

The first day started with a talk by Jörg Kienitz. He described the recent research on SABR models and the ways that people are making them work with negative rates (one of the several themes of the meeting). There was a lot of stuff covered; you'll find it in the slides when they're available. As Jörg mentioned, QuantLib contains an implementation of the SABR and ZABR models, even though it's not up to date with all the recent developments.

The second talk was from Alexander Sokol, describing the work done by his company on the tapescript library. It provides a numerical type meant to be a drop-in replacement for double in AAD calculations; as a proof of concept, they managed to recompile the whole QuantLib using it. It's an all-or-nothing approach that can be used to quickly convert an existing legacy library to AAD. The typescript library is available on GitHub.

After lunch, the afternoon started with a joint talk by Ferdinando Ametrano and Paolo Mazzocchi describing their results on modeling tenor basis spread between overnight and forward curves. As it turns out, an abcd parameterization provides quite a nice fit of the observed basis and can be made an exact fit with some small correction coefficients. The approach looks promising, the slides are already published, and an initial implementation is available on Paolo's fork of QuantLib on GitHub (look for the new code in the ql/experimental/tenorbasis folder). A paper will follow.

Next, a more technological talk. Eric Ehlers reported on the status of his reposit project, the successor to ObjectHandler to be used in future versions of QuantLibXL. It looks good. Also, Eric left the podium for part of his presentation to Cristian Alzati from Sayula, who demonstrated their joint work on the =countify platform, making Excel spreadsheets available on the cloud. QuantLibXL is one of the available addins.

In the final talk of the day, Klaus Spanderen and Johannes Göttker-Schnetmann presented the conclusion of the work, previewed in last year's meeting, on the calibration of stochastic local volatilities. Their code is available as a pull request, and will be probably merged in the next QuantLib release.

And then there was beer and a well-deserved rest.

The second day was opened by Peter Caspers, who made two short presentations because that's how he rolls. The first was on AAD, and described his efforts in enabling part of QuantLib to use it. His approach complements that of Alexander, and trades ease of conversion for better performance. You can read about it on his blog, too. I'll be watching closely for a possible synergy of the two approaches; in any case, it's nice to have the choice. In his second presentation, Peter reviewed the work done to enable QuantLib work with negative rates. In short: we're almost there, but there's still a couple of pull requests missing.

Andreas Pfadler, the next speaker, had the best opening line of the meeting. Quoting from memory, it was something like "this work is the fruit of many lonely nights in a hotel". Andreas reported on his development of an open-source architecture for distributing calculations, using the QuantLib Java bindings among other tools and adding scripting capabilities via Scala. During his demo, he also overcame a technical glitch on his computer. Another speaker later confessed over coffee that he would have been crushed by it. Andreas took it in stride by moving the demo to the Amazon cloud.

In the last talk of the morning, Roland Lichters fought bravely with a sore throat to tell us about CSA pricing using QuantLib, especially in the context of the added complexity and the pricing changes caused by using negative rates for collateral discounting. The details are on the slides, and caused an interesting discussion afterward since there's not a lot of agreement on how the effects should be modeled.

In the afternoon, another couple of talks. In the first, Sebastian Schlenkrich came back to the topic of tenor basis spreads; his take was on transforming rate volatilities in this kind of models. Again, you'll be able to read all about it in his slides when they're available. In the second one, I described shortly how I used the IPython Notebook and Docker together with QuantLib. Stay tuned on this channel for more information.

And with that, we were off to our flights. Again, thanks to the organizers, the sponsors, the speakers, and all the participants. Here's to meeting again next year.

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

## Thursday, November 26, 2015

### A quick look at the QuantLib 1.7 release

Hello everybody.

Last Monday, I released QuantLib 1.7 (you can download it here if you still haven't). As is becoming a tradition, here's a quick look at the release.

It comes five months after 1.6, which is not bad. (In the meantime, we also had two small releases in the 1.6 series for compatibility with the latest Boost and Visual C++.) It contains 53 issues and pull requests that you can examine in detail on GitHub; just follow the previous link. The usual git incantation shows that the release consists of 217 commits by 14 people (once you filter out a few duplicates).

Other people were involved by contributing bug reports or patches; their names, as far as I could find them, are in the list of change for the release. I hope I didn't miss anyone; if so, I apologize.

There are two new features I'd like to point out in particular. Both are disabled by default, as they cause a performance penalty. The first is the addition of a time to the Date class; this makes it possible to finally price intraday options. The second is a reimplementation of the Observer pattern that makes it safe to use from C# or Java, whose garbage collector had a habit of sometime destroying objects during notification and crashing the whole thing down. Both new features are mostly the work of Klaus Spanderen. You're encouraged to try them out; the release notes explain how to enable them (look for them at the bottom of the file).

That's all for this post. See you after Düsseldorf.

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

## Thursday, November 19, 2015

### Chapter 8, part 6 of n: example, American option

Welcome back.

This week, I finally continue with content from chapter 8 (the one on the finite-difference framework, a.k.a. the one that's taking forever).

Next post will probably be a report from the QuantLib User Meeting, in two or three weeks. There's probably still places available, so click here for more info.

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

### Example: American option

At this point, writing a finite-difference pricing engine should be just a matter of connecting the dots. Well, not quite. In this section, I'll sketch the implementation of an American-option engine in QuantLib, which is somewhat more complex than expected (as you can see for yourself from the following figure.)

My reasons for doing this are twofold. On the one hand, I nearly got lost myself when I set out to write this chapter and went to read the code of the engine; so I thought it might be useful to put a map out here. On the other hand, this example will help me draw a comparison with the new (and more modular) framework.

Mind you, I'm not dissing the old implementation. The reason it got so complex was that we tried to abstract out reusable chunks of code, which makes perfect sense. The problem is that, although we didn't see it at the time, inheritance was probably the wrong way to do it.

Let's start with the FDVanillaEngine class, shown in the listing below. It can be used as a base class for both vanilla-option and dividend vanilla-option engines, which might explain why the name is not as specific as, say, FDVanillaOptionEngine. (We might just have decided to shorten the name, though. I don't think anybody remembers after all these years.)
    class FDVanillaEngine {
public:
FDVanillaEngine(
const shared_ptr<GeneralizedBlackScholesProcess>&,
Size timeSteps, Size gridPoints,
bool timeDependent = false);
virtual ~FDVanillaEngine() {}
protected:
virtual void setupArguments(
const PricingEngine::arguments*) const;
virtual void setGridLimits() const;
virtual void initializeInitialCondition() const;
virtual void initializeBoundaryConditions() const;
virtual void initializeOperator() const;

shared_ptr<GeneralizedBlackScholesProcess> process_;
Size timeSteps_, gridPoints_;
bool timeDependent_;
mutable Date exerciseDate_;
mutable boost::shared_ptr<Payoff> payoff_;
mutable TridiagonalOperator finiteDifferenceOperator_;
mutable SampledCurve intrinsicValues_;
typedef BoundaryCondition<TridiagonalOperator> bc_type;
mutable std::vector<boost::shared_ptr<bc_type> > BCs_;

virtual void setGridLimits(Real, Time) const;
virtual Time getResidualTime() const;
void ensureStrikeInGrid() const;
private:
Size safeGridPoints(Size gridPoints,
Time residualTime) const;
};

This class builds most of the pieces required for a finite-difference model, based on the data passed to its constructor: a Black-Scholes process for the underlying, the number of desired time steps and grid points, and a flag that I'm going to ignore until the next subsection. Besides the passed inputs, the data members of the class include information to be retrieved from the instrument (that is, the exercise date and payoff) and the pieces of the model to be built: the differential operator, the boundary conditions, and the array of initial values corresponding to the intrinsic values of the payoff. The latter array is stored in an instance of the SampledCurve class, which adds a few utility methods to the stored data.

The rest of the class interface is made of protected methods that builds and operates on the data members. I'll just go over them quickly: you can read their implementation in the library for more details.

First, the spectacularly misnamed setupArguments method does the opposite of its namesake in the Instrument class: it reads the required exercise and payoff information from the passed arguments structure and copies them into the corresponding data members of FDVanillaEngine.

The setGridLimits method determines and stores the minimum and maximum value of the logarithmic model grid, based on the variance of the passed process over the residual time of the option. The calculation enforces that the current value of the underlying is at the center of the grid, that the strike value is within its range, and that the number of its points is large enough. (I'd note that the method might override the number of grid points passed by the user. In hindsight, I'm not sure that doing it silently is a good idea.) The actual work is delegated to a number of other methods: an overloaded version of setGridLimits, safeGridPoints, and ensureStrikeInGrid.

The initializeInitialCondition method fills the array of intrinsic values by sampling the payoff on the newly specified grid; thus, it must be called after the setGridLimits method.

The initializeBoundaryConditions method, to be called as the next step, instantiates the lower and upper boundary conditions. They're both Neumann conditions, and the value of the derivative to be enforced is calculated numerically from the array of intrinsic values.

Finally, the initializeOperator method creates the tridiagonal operator based on the calculated grid and the stored process. Again, the actual work is delegated: in this case, to the OperatorFactory class, that I'll describe later.

All of these methods are declared as virtual, so that the default implementations can be overridden if needed. This is not optimal: in order to change any part of the logic one has to use inheritance, which introduces an extra concept just for customization and doesn't lend itself to different combinations of changes. A Strategy pattern would be better, and would also make some of the logic more reusable by other instruments.

All in all, though, the thing is manageable: see the FDEuropeanEngine class, shown in the next listing, which can be implemented in a reasonable amount of code.
    template <template <class> class Scheme = CrankNicolson>
class FDEuropeanEngine : public OneAssetOption::engine,
public FDVanillaEngine {
public:
FDEuropeanEngine(
const shared_ptr<GeneralizedBlackScholesProcess>&,
Size timeSteps=100, Size gridPoints=100,
bool timeDependent = false);
private:
mutable SampledCurve prices_;
void calculate() const {
setupArguments(&arguments_);
setGridLimits();
initializeInitialCondition();
initializeOperator();
initializeBoundaryConditions();

FiniteDifferenceModel<Scheme<TridiagonalOperator> >
model(finiteDifferenceOperator_, BCs_);

prices_ = intrinsicValues_;

model.rollback(prices_.values(), getResidualTime(),
0, timeSteps_);

results_.value = prices_.valueAtCenter();
results_.delta = prices_.firstDerivativeAtCenter();
results_.gamma = prices_.secondDerivativeAtCenter();
results_.theta = blackScholesTheta(process_,
results_.value,
results_.delta,
results_.gamma);
}
};

Its calculate method sets everything up by calling the appropriate methods from FDVanillaEngine, creates the model, starts from the intrinsic value of the option at maturity and rolls it back to the evaluation date. The value and a couple of Greeks are extracted by the corresponding methods of the SampledCurve class, and the theta is calculated from the relationship that the Black-Scholes equation imposes between it and the other results. (By replacing the derivatives with the corresponding Greeks, the Black-Scholes equation says that $$\Theta + \frac{1}{2} \sigma^2 S^2 \Gamma + (r-q)S\Delta -rV = 0$$.)

What with the figure above then? How do we get three levels of inheritance between FDVanillaEngine and FDAmericanEngine? It was due to the desire to reuse whatever pieces of logic we could. As I said, the idea was correct: there are other options in the library that use part of this code, such as shout options, or options with discrete dividends. The architecture could be improved.

First, we have the FDStepConditionEngine, sketched in the following listing.
    template <template <class> class Scheme = CrankNicolson>
class FDStepConditionEngine : public FDVanillaEngine {
public:
FDStepConditionEngine(
const shared_ptr<GeneralizedBlackScholesProcess>&,
Size timeSteps, Size gridPoints,
bool timeDependent = false);
protected:
// ...data members...
virtual void initializeStepCondition() const = 0;
virtual void calculate(PricingEngine::results*) const {
OneAssetOption::results * results =
dynamic_cast<OneAssetOption::results *>(r);
setGridLimits();
initializeInitialCondition();
initializeOperator();
initializeBoundaryConditions();
initializeStepCondition();

typedef /* ... */ model_type;

prices_ = intrinsicValues_;
controlPrices_ = intrinsicValues_;
// ...more setup (operator, BC) for control...

model_type model(operatorSet, bcSet);
model.rollback(arraySet, getResidualTime(),
0.0, timeSteps_, conditionSet);

results->value = prices_.valueAtCenter()
- controlPrices_.valueAtCenter()
+ black.value();
// same for Greeks
}
};

It represent a finite-difference engine in which a step condition is applied at each step of the calculation. In its calculate method, it implements the bulk of the pricing logic—and then some. First, it sets up the data members by calling the methods inherited from FDVanillaEngine, as well as an initializeStepCondition method that it declares as pure virtual and that derived classes must implement: it must create an instance of the StepCondition class appropriate for the given engine. Then, it creates two arrays of values; the first for the option being priced, and the second for a European option that will be used as a control variate (this also requires to set up a corresponding operator, as well as a pricer object implementing the analytic Black formula). Finally, the model is created and used for both arrays, with the step condition being applied only to the first one, and the results are extracted and corrected for the control variate.

I don't have particular bones to pick with this class, except for the name, which is far too generic. I'll just add a note on the usage of control variates. We have already seen the technique in a previous post, where it was used to narrow down the width of the simulated price distribution; here it is used to improve the numerical accuracy. It is currently forced upon the user, since there's no flag allowing to enable or disable it; and it is relatively more costly than in Monte Carlo simulations (there, the path generation is the bulk of the computation and is shared between the option and the control; here, using it almost doubles the computational effort). The decision of whether it's worth using should be probably be left to the user. Also, we should use temporary variables for the control data instead of declaring other mutable data members; they're turning into a bad habit.

Next, the FDAmericanCondition class template, shown in the next listing.
    template <typename baseEngine>
class FDAmericanCondition : public baseEngine {
public:
FDAmericanCondition(
const shared_ptr<GeneralizedBlackScholesProcess>&,
Size timeSteps = 100, Size gridPoints = 100,
bool timeDependent = false);
protected:
void initializeStepCondition() const;
};

It takes its base class as a template argument (in our case, it will be FDVanillaEngine) and provides the initializeStepCondition method, which returns an instance of the AmericanCondition class. Unfortunately, the name FDAmericanCondition is quite confusing: it suggests that the class is a step condition, rather than a building block for a pricing engine.

The next to last step is the FDEngineAdapter class template.
    template <typename base, typename engine>
class FDEngineAdapter : public base, public engine {
public:
const shared_ptr<GeneralizedBlackScholesProcess>& p,
Size timeSteps=100, Size gridPoints=100,
bool timeDependent = false)
: base(p, timeSteps, gridPoints,timeDependent) {
this->registerWith(p);
}
private:
void calculate() const {
base::setupArguments(&(this->arguments_));
base::calculate(&(this->results_));
}
};

It connects an implementation and an interface by taking them as template arguments and inheriting from both: in this case, we'll have FDAmericanCondition as the implementation and OneAssetOption::engine as the interface. The class also provides a bit of glue code in its calculate method that satisfies the requirements of the engine interface by calling the methods of the implementation.

Finally, the FDAmericanEngine class just inherits from FDEngineAdapter and specifies the classes to be used as bases.
    template <template <class> class Scheme = CrankNicolson>
class FDAmericanEngine
FDAmericanCondition<
FDStepConditionEngine<Scheme> >,
OneAssetOption::engine> {
public:
FDAmericanEngine(
const shared_ptr<GeneralizedBlackScholesProcess>&,
Size timeSteps=100, Size gridPoints=100,
bool timeDependent = false);
};

The question is whether it is worth to increase the complexity of the hierarchy in order to reuse the bits of logic in the base classes. I'm not sure I have an answer, but I can show an alternate implementation and let you make the comparison on your own. If we let FDAmericanEngine inherit directly from FDStepConditionEngine and OneAssetOption::engine, and if we move into this class the code from both FDAmericanCondition and FDEngineAdapter (that we can remove afterwards), we obtain the implementation in the listing below.
    template <template <class> class Scheme = CrankNicolson>
class FDAmericanEngine
: public FDStepConditionEngine<Scheme>,
public OneAssetOption::engine {
typedef FDStepConditionEngine<Scheme> fd_engine;
public:
FDAmericanEngine(
const shared_ptr<GeneralizedBlackScholesProcess>& p,
Size timeSteps=100, Size gridPoints=100,
bool timeDependent = false)
: fd_engine(p, timeSteps, gridPoints, timeDependent) {
this->registerWith(p);
}
protected:
void initializeStepCondition() const;
void calculate() const {
fd_engine::setupArguments(&(this->arguments_));
fd_engine::calculate(&(this->results_));
}
};

My personal opinion? I tend to lean towards simplicity in my old age. The code to be replicated would be little, and the number of classes that reuse it is not large (about half a dozen in the current version of the library). Moreover, the classes that we'd remove (FDAmericanCondition and FDEngineAdapter) don't really model a concept in the domain, so I'd let them go without any qualms. Too much reuse without a proper abstraction might be a thing, after all.

A final note: as you can see, in this framework there are no high-level classes encapsulating generic model behavior, such as McSimulation for Monte Carlo (see this post). Whatever logic we had here was written in classes meant for a specific instrument—in this case, plain options in a Black-Scholes-Merton model

## Thursday, November 12, 2015

### Upcoming: QuantLib User Meeting 2015

Hello everybody.

It's been a while, I know. October has been a busy month.

Just a quick post to remind you of a bit of news you might have read on the QuantLib mailing list: IKB and CompatibL are teaming up to sponsor the next QuantLib user meeting in Düsseldorf. Like the past two meetings (which were excellent: see my posts on them here and here) it will span two days, namely, November 30th and December 1st. Attendance is free but places are limited, so hurry up if you're interested (you should). More info at this link.

And that's all, folks. I'll see you soon with some new content.

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

## Tuesday, September 29, 2015

### QuantLib notebook: term structures and reference dates

Hello everybody.

Well, it's been a while. Last time I posted, we were still in August (and there were a number of news in that post, so you might want to go back and read it if you missed it). I have since released QuantLib 1.6.2, which you'll need if you want to use Visual Studio and Boost 1.59 and you can disregard otherwise. And hey, look at that: the Google+ buttons at the bottom of the post have a new look, too. What will they think of next?

This week, a new QuantLib notebook: number 6 in an ongoing series, if you're keeping count. These screencast could turn into a book, so go to this Leanpub page and tell them you're interested if you want this to happen (and while you're there, have a look at Implementing QuantLib, too).

Also, let me bother you for half a minute more: there are still places left for my next Introduction to QuantLib Development course. Click through for details.

Ok, here's the notebook. Set it to full screen, sit back and enjoy.

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

## Monday, August 24, 2015

### Chapter 8, part 5 of n: finite-difference models

Welcome back, everybody.

This week, the next part of the series on the finite-difference framework.

But first, a small personal milestone: while I was in vacation, the 100th copy of Implementing QuantLib was sold on Leanpub. So, to the 100th reader who bought it between last Friday and Saturday, but also to all the others: thank you.

A couple of weeks ago we also had... well, actually, I don't know if we had some kind of event or not. One morning I found in my inbox a few hundreds notifications that email addresses were removed from the QuantLib mailing lists. This usually happens when an address bounces a few times; the mailing-list software then flags it as obsolete and removes it. However, we never had that many at once. I'm not sure if this is some kind of glitch after the recent Sourceforge outages, or the result of their changing some kind of threshold. In any case, if you haven't received emails from the list in a while, please check that you're still subscribed. You can do so from http://quantlib.org/mailinglists.shtml: follow the "Subscribe/Unsubscribe/Preferences" link for the lists you're subscribed to.

In the "back to business" department: there's a few things happening in October. If you're going to the WBS Fixed Income Conference in Paris, you may find me there. On Wednesday 7trh, I'll be in the bench during Alexander Sokol's workshop on how to add AAD to QuantLib (I'm listed as a presenter, but I'll be just there ready to intervene if bits of particular QuantLib knowledge are needed). On the morning of Friday 9th, instead, I'll have a slot to present QuantLib and its next developments in fixed income.

Also in October, I'll be running my next Introduction to QuantLib Development course (in London, from Monday 19th to Wednesday 21st). Details and a registration form are at http://bit.ly/quantlib2015; you can also check out the hangout I did with Jacob Bettany to talk about it.

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

### The FiniteDifferenceModel class

The FiniteDifferenceModel class, shown in the next listing, uses an evolution scheme to roll an asset back from a later time $$t_1$$ to an earlier time $$t_2$$, taking into account any conditions we add to the model and stopping at any time we declare as mandatory (because something is supposed to happen at that point).

A word of warning: as I write this, I'm looking at this class after quite a few years and I realize that some things might have been done differently—not necessarily better, mind you, but there might be alternatives. After all this time, I can only guess at the reasons for such choices. While I describe the class, I'll point out the alternatives and write my guesses.
    template<class Evolver>
class FiniteDifferenceModel {
public:
typedef typename Evolver::traits traits;
typedef typename traits::operator_type operator_type;
typedef typename traits::array_type array_type;
typedef typename traits::bc_set bc_set;
typedef typename traits::condition_type condition_type;

FiniteDifferenceModel(
const Evolver& evolver,
const std::vector<Time>& stoppingTimes =
std::vector<Time>());

FiniteDifferenceModel(
const operator_type& L,
const bc_set& bcs,
const std::vector<Time>& stoppingTimes =
std::vector<Time>());
const Evolver& evolver() const{ return evolver_; }
void rollback(array_type& a,
Time from,
Time to,
Size steps) {
rollbackImpl(a,from,to,steps,
(const condition_type*) 0);
}
void rollback(array_type& a,
Time from,
Time to,
Size steps,
const condition_type& condition) {
rollbackImpl(a,from,to,steps,&condition);
}
private:
void rollbackImpl(array_type& a,
Time from,
Time to,
Size steps,
const condition_type* condition);
Evolver evolver_;
std::vector<Time> stoppingTimes_;
};

I'll just mention quickly that the class takes the evolution scheme as a template argument (as in the code, I might call it evolver for brevity from now on) and extracts some types from it via traits; that's a technique you saw a few times already in earlier chapters.

So, first of all: the constructors. As you see, there's two of them. One, as you'd expect, takes an instance of the evolution scheme and a vector of stopping times (and, not surprisingly, stores copies of them). The second, instead, takes the raw ingredients and builds the evolver.

The reason for the second one might have been convenience, that is, to avoid writing the slighly redundant
    FiniteDifferenceModel<SomeScheme>(SomeScheme(L, bcs), ts);

but I'm not sure that the convenience offsets having more code to maintain (which, by the way, only works for some evolvers and not for others; for instance, the base MixedScheme class doesn't have a constructor taking two arguments). This constructor might make more sense when using one of the typedefs provided by the library, such as StandardFiniteDifferenceModel, which hides the particular scheme used; but I'd argue that you should know what scheme you're using anyway, since they have different requirements in terms of grid spacing and time steps. (By the same argument, the StandardFiniteDifferenceModel typedef would be a suggestion of what scheme to use instead of a way to hide it.)

The main logic of the class is coded in the rollback and rollbackImpl methods. The rollback method is the public interface and has two overloaded signatures: one takes the array of values to roll back, the initial and final times, and the number of steps to use, while the other also takes a StepCondition instance to use at each step. Both forward the call to the rollbackImpl method, which takes a raw pointer to step condition as its last argument; the pointer is null in the first case and contains the address of the condition in the second.

Now, this is somewhat unusual in the library: in most other cases, we'd just have a single public interface taking a shared_ptr and defining a null default value. I'm guessing that this arrangement was made to allow client code to simply create the step condition on the stack and pass it to the method, instead of going through the motions of creating a shared pointer (which, in fact, only makes sense if it's going to be shared; not passed, used and forgotten). Oh, and by the way, the use of a raw pointer in the private interface doesn't contradict the ubiquitous use of smart pointers everywhere in the library. We're not allocating a raw pointer here; we're just passing the address of an object living on the stack (or inside a smart pointer, for that matters) and which can manage its lifetime very well, thank you.

One last remark on the interface: we're treating boundary conditions and step conditions differently, since the former are passed to the constructor of the model and the latter are passed to the rollback method. (Another difference is that we're using explicit sets of boundary conditions but just one step condition, which can be a composite of several ones if needed. I have no idea why.) One reason could have been that this allows more easily to use different step conditions in different calls to rollback, that is, at different times; but that could be done inside the step condition itself, as it is passed the current time. The real generalization would have been to pass the boundary conditions to rollback, too, as it's currently not possible to change them at a given time in the calculation. This would have meant changing the interfaces of the evolvers, too: the boundary conditions should have been passed to their step method, instead of their constructors.

And finally, the implementation.
    void FiniteDifferenceModel::rollbackImpl(
array_type& a,
Time from,
Time to,
Size steps,
const condition_type* condition) {
Time dt = (from-to)/steps, t = from;
evolver_.setStep(dt);
if (!stoppingTimes_.empty() && stoppingTimes_.back()==from) {
if (condition)
condition->applyTo(a,from);
}
for (Size i=0; i<steps; ++i, t -= dt) {
Time now = t, next = t-dt;
if (std::fabs(to-next) < std::sqrt(QL_EPSILON))
next = to;
bool hit = false;
for (Integer j = stoppingTimes_.size()-1; j >= 0 ; --j) {
if (next <= stoppingTimes_[j]
&& stoppingTimes_[j] < now) {
hit = true;
evolver_.setStep(now-stoppingTimes_[j]);
evolver_.step(a,now);
if (condition)
condition->applyTo(a,stoppingTimes_[j]);
now = stoppingTimes_[j];
}
}
if (hit) {
if (now > next) {
evolver_.setStep(now - next);
evolver_.step(a,now);
if (condition)
condition->applyTo(a,next);
}
evolver_.setStep(dt);
} else {
evolver_.step(a,now);
if (condition)
condition->applyTo(a, next);
}
}
}

We're going back from a later time from to an earlier time to in a given number of steps. This gives us a length dt for each step, which is set as the time step for the evolver. If from itself is a stopping time, we also apply the step condition (if one was passed; this goes for all further applications). Then we start going back step by step, each time going from the current time t to the next time, t-dt; we'll also take care that the last time is actually to and possibly snap it back to the correct value, because floating-point calculations might have got us almost there but not quite.

If there were no stopping times to hit, the implementation would be quite simple; in fact, the whole body of the loop would be the three lines inside the very last else clause:
    evolver_.step(a,now);
if (condition)
condition->applyTo(a, next);

that is, we'd tell the evolver to roll back the values one step from now, and we'd possibly apply the condition at the target time next. This is what actually happens when there are no stopping times between now and next, so that the control variable hit remains false.

In case there are one or more stopping times in the middle of the step, things get a bit more complicated. Have a look at the upper half of the following figure, where I sketched a series of regular steps with two stopping times, $$s_1$$ and $$s_2$$, falling in the middle of two of them.

Let's say we have just made the step from $$t_7$$ to $$t_6$$ as described in the previous paragraph, and should now go back to $$t_5$$. However, the code will detect that it should stop at $$s_2$$, which is between now and next; therefore, it will set the control variable hit to true, it will set the step of the evolver to the smaller step $$t_6-s_2$$ that will bring it to the stopping time, and will perform the step (again, applying the condition afterwards). Another, similar change of step would happen if there were two or more stopping times in a single step; the code would then step from one to the other. Finally, the code will enter the if (hit) clause, which get things back to normal: it performs the remaining step $$s_2-t_5$$ and then resets the step of the evolver to the default value, so that it's ready for the next regular step to $$t_4$$.

As a final remark, let's go back to trees for a minute. Talking about the tree framework, I mentioned in passing that the time grid of a lattice is built to include the mandatory stopping times. In this case, the grid would have been the one shown in the lower half of the previous figure. As you see, the stopping times $$s_1$$ and $$s_2$$ were made part of the grid by using slightly smaller steps between to and $$s_1$$ and slightly larger steps between $$s_2$$ and from. (In this case, the change results in the same number of steps as the finite-difference grid. This is not always guaranteed; the tree grid might end up having a couple of steps more or less than the requested number.) I'm not aware of drawbacks in using either of the two methods, but experts might disagree so I'm ready to stand corrected. In any case, should you want to replicate the tree grid in a finite-difference model, you can do so by making several calls to rollback and going explicitly from stopping time to stopping time.

## Monday, July 27, 2015

### Chapter 8, part 4 of n: step conditions

Welcome back.

The content for this post is the fourth part of an ongoing series on the QuantLib finite-difference framework.

But first: this week the final version of Visual Studio 2015 was released and, whaddayaknow, QuantLib 1.6 was not quite prepared for it. We had tried to be forward-looking and had prepared our version-detecting code for VC13, but alas, Visual Studio 2015 turned out to be VC14 instead. The library compiles and links just fine, but has the wrong name and the auto-link pragmas get confused. This week I might be releasing a 1.6.1 version that will fix the issue, God willing and the creek don't rise—since the other bit of news this week was that Sourceforge had a major outage (you might have noticed that the QuantLib web site and the mailing lists have been down for a few days). Now it's mostly up again, but as of today it's still not possible to make new releases. Keep your fingers crossed.

Oh, and my Quants Hub workshop is still a steal at £99, but from what I'm told this will only last for a few more days, until the end of July. Get it cheap while you can.

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

### Step conditions

At any step (or in the middle of a step; we'll see later how we account for this) something might happen that changes the price of the asset; for instance, a right can be exercised or a coupon could be paid.

As opposed to boundary conditions, we called these step conditions. The name might not be the best; I can easily see the model as enforcing the "condition" that the asset price is the maximum between its continuation and the intrinsic value, but a coupon payment is harder to categorize in this way. In any case, the base class for all such conditions is shown in the listing below.
    template <class array_type>
class StepCondition {
public:
virtual ~StepCondition() {}
virtual void applyTo(array_type& a, Time t) const = 0;
};

The applyTo method, which modifies the array of asset values in place, must also take care of checking that the condition applies at the given time; so, for instance, an American exercise will change the asset values regardless of the time, whereas a dividend payment will first check the passed time and bail out if it doesn't match the known payment time.

As you see, the base class is quite simple; but derived classes over-generalized and went South. In the next listing, you can see one such class, from which even simple conditions (like the one for American exercise) are derived. We could have implemented a simple applyTo method comparing the option values and the intrinsic values. Instead, we started adding constructors over the years so that instances could be built from either an existing array of intrinsic values, or a payoff, or some data that could describe the payoff. Of course these alternatives needed to be stored into different data members, so we added type erasure to the mix, we defined an interface to retrieve intrinsic values, and we gave it an array-based and a payoff-based implementation.
    template
class CurveDependentStepCondition
: public StepCondition {
public:
void applyTo(Array &a, Time) const {
for (Size i = 0; i < a.size(); i++) {
a[i] = applyToValue(a[i], getValue(a,i));
}
}
protected:
CurveDependentStepCondition(Option::Type type,
Real strike);
CurveDependentStepCondition(const Payoff *p);
CurveDependentStepCondition(const array_type & a);

class CurveWrapper {
public:
Real getValue(const array_type &a,
Size index) const = 0;
};
boost::shared_ptr curveItem_;

Real getValue(const array_type &a, Size index) const {
return curveItem_->getValue(a, index);
}

virtual Real applyToValue(Real, Real) const = 0;

class ArrayWrapper : public CurveWrapper {
array_type value_;
public:
...
};

class PayoffWrapper : public CurveWrapper {
boost::shared_ptr payoff_;
public:
...
};
};

Ironically, we only ever used half of this. Existing engines use the array-based constructor and implementation, and the payoff-based inner class is not used anywhere in the library, Which is just as well; because, as I reviewed the code to write this chapter, I found out that its implementation is wrong. So all in all, we could have written a much simpler implementation of an American exercise condition that would just take an array of intrinsic values and implement applyTo accordingly. If you ever implement a step condition, I suggest you do the same.

A final note; the step condition is not applied inside the step method of the evolution scheme, as we'll see shortly. What this means is that there's no guarantee that the step condition won't break the boundary condition, or even the other way around. In practice, though, applying a step condition for a model that makes any sense will naturally enforce the chosen boundary condition.

Next time, with all the basic pieces in place, we'll start to put them together.

#### Aside: look, Ma, no hands

The faulty implementation of PayoffWrapper was added in February 2003, which probably makes this the longest-lived bug in the library (apart from those we still haven't found, of course). If I hadn't been looking at the code to write this chapter, the critter could have got old enough to drive.

Now, bugs happen; nothing shocking about this. What irks me, though, is that the problem would have been caught easily by a couple of unit tests. I and others usually try and review what goes in the library, and the code that gets contributed has often been used in the field and purged of obvious errors; but still I'm somewhat worried about the lack of coverage of old code in the test suite. If you have a coverage tool and some time on your hands, please let me hear from you. I'd love to see some analysis on this.

## Monday, July 20, 2015

### MoneyScience Hangout available on YouTube

Hello everybody.

A short bit of news: that Google Hangout I told you I would do with MoneyScience's Jacob Bettany? We did it, and it's now online on YouTube. Here it is. More details on my next course (which will be in London on October 19th to 21st) are at http://bit.ly/quantlib2015.

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

## Monday, July 13, 2015

### Chapter 8, part 3 of n: boundary conditions

Hello everybody.

Just a reminder before this week's content (that is, another installment in the series on finite-difference methods): my Quants Hub workshop, A Look at QuantLib Usage and Development, is going to be £99 until the end of July. Grab it cheap while you can.

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

### Boundary conditions

As I said, there's more going on in the step method. To begin with, we'll have to enforce boundary conditions at either end of the grid; for instance, we might want to ensure that an option has price close to 0 at the grid boundary which is deep out of the money and derivative close to 1 at the boundary which is deep in the money. Such conditions are modeled by the BoundaryCondition class, shown in the listing below.
    template <class Operator>
class BoundaryCondition {
public:
typedef Operator operator_type;
typedef typename Operator::array_type array_type;

enum Side { None, Upper, Lower };

virtual ~BoundaryCondition() {}

virtual void applyBeforeApplying(operator_type&) const = 0;
virtual void applyAfterApplying(array_type&) const = 0;
virtual void applyBeforeSolving(operator_type&,
array_type& rhs) const = 0;
virtual void applyAfterSolving(array_type&) const = 0;
};

As you see, the class uses runtime polymorphism. This is because we wanted to store collections of them, possibly of different types (you might remember the bc_set typedef from the listing of the MixedScheme class), and most suitable containers are homogeneous and thus require a common base class (std::tuple wasn't even a gleam in the standard committee's eye, and handling it requires template metaprogramming, which I'd try to avoid unless necessary—even in this time and age) We could have used std::pair if we had decided to limit ourselves to the 1-D case; we didn't, though, and in hindsight I think it avoided adding more complexity to the framework.

The interface of the class is, well, peculiar. Glossing over the Side enumeration (which is meant to specify a grid side in 1-D models, and thus too specific), methods like applyBeforeApplying are plausible contenders for the worst name in the library. The idea here is that the boundary condition can be enforced, or "applied", in different ways. One can wait for the operator's applyTo or solveFor to execute, and then modify the result; this would be implemented by applyAfterApplying and applyAfterSolving, respectively. Another possibility is to set up the operator and the input array beforehand, so that the result of the operation will satisfy the condition; this is done in applyBeforeApplying and applyBeforeSolving. (For some reason, the former only takes the operator. It might have been an oversight, due to the fact that the existing conditions didn't require to modify the array.)

As an example, look at the DirichletBC class in the listing below. It implements a simple Dirichlet boundary condition, i.e., one in which the function value at a given end of the grid must equal a given constant value.
    class DirichletBC
: public BoundaryCondition<TridiagonalOperator> {
public:
DirichletBC(Real value, Side side);

void applyBeforeApplying(TridiagonalOperator& L) const {
switch (side_) {
case Lower:
L.setFirstRow(1.0,0.0);
break;
case Upper:
L.setLastRow(0.0,1.0);
break;
}
}
void applyAfterApplying(Array& u) const {
switch (side_) {
case Lower:
u[0] = value_;
break;
case Upper:
u[u.size()-1] = value_;
break;
}
}
void applyBeforeSolving(TridiagonalOperator& L,
Array& rhs) const {
switch (side_) {
case Lower:
L.setFirstRow(1.0,0.0);
rhs[0] = value_;
break;
case Upper:
L.setLastRow(0.0,1.0);
rhs[rhs.size()-1] = value_;
break;
}
}
void applyAfterSolving(Array&) const {}
};

As you can see, it is no longer a class template; when implemented, boundary conditions are specialized for a given operator—not surprisingly, because they have to know how to access and modify it. In this case, we're using the tridiagonal operator I described in an earlier subsection.

The constructor is simple enough: it takes the constant value and the side of the grid to which it must be applied, and stores them.

As I said, the applyBeforeApplying method must set up the operator L so that the result of L.apply(u), or $$\mathbf{u^{\prime}} = L \cdot \mathbf{u}$$, satisfies the boundary condition. To do this, it sets the first (last) row to the corresponding row of the identity matrix; this ensures that the first (last) element of $$\mathbf{u^{\prime}}$$ equal the corresponding element of $$\mathbf{u}$$. Given that the array satisfied the boundary condition at the previous step, it keeps satisfying it at this step. Relying on this is a bit unsafe (for instance, it would break if the value changed at different times, or if any other event modified the array values) but it's the best one can do without accessing the input array. The applyAfterApplying method is simpler and safer: it just sets the value of the output array directly.

The applyBeforeSolving method works as its sibling applyBeforeApplying, but it can also access the input array, so it can ensure that it contains the correct value. Finally, and surprisingly enough, as of this writing the applyAfterSolving method is empty. When we wrote this class, we probably relied on the fact that applyBeforeSolving would also be called, which doesn't strike me as a very good idea in hindsight: we're not sure of what methods any given evolution scheme may call (this also makes applyBeforeApplying less safe than it seems).

The reason this worked is the way the MixedScheme::step method is implemented, as seen in the next listing. At each step, all boundary conditions are enforced first before and after applying the explicit part of the scheme, and then before and after solving for the implicit part. In the case of the Dirichlet condition, this causes applyAfterApplying to fix any problem that applyBeforeApplying might have introduced and also ensures that applyBeforeSolving does the work even if applyAfterSolving doesn't.

If a boundary condition class implemented all of its methods correctly, this would be redundant; enforcing the condition just after the operations would be enough. However, we probably assumed that some conditions couldn't (or just didn't) implement all of them, and went for the safe option.
    template <class Operator>
void MixedScheme<Operator>::step(array_type& a, Time t) {
if (theta_ != 1.0) { // there is an explicit part
for (Size i=0; i<bcs_.size(); i++)
bcs_[i]->applyBeforeApplying(explicitPart_);
a = explicitPart_.applyTo(a);
for (Size i=0; i<bcs_.size(); i++)
bcs_[i]->applyAfterApplying(a);
}
if (theta_ != 0.0) { // there is an implicit part
for (Size i=0; i<bcs_.size(); i++)
bcs_[i]->applyBeforeSolving(implicitPart_, a);
implicitPart_.solveFor(a, a);
for (Size i=0; i<bcs_.size(); i++)
bcs_[i]->applyAfterSolving(a);
}
}

If we ever set out to revamp this part of the library (not that likely, given the new framework) this is something we should get a look at. The problem to solve, unfortunately, is not an easy one: to call only the methods that are required by the scheme and supported by the boundary condition, we'd need to know the specific type of both.

Inverting the dependency and passing the scheme to the boundary condition, instead of the other way around, wouldn't work: the boundary condition wouldn't know if the scheme is calling applyTo or solveFor, and wouldn't be able to act inside the scheme's step method (as in MixedScheme::step, where the boundary conditions are enforced between the implicit and explicit parts of the step.

A mediator or some kind of visitor might work, but the complexity of the code would increase (especially if the number of conditions grow). Had they been available in C++, multimethods might have been the best choice; but even those might not be usable for multiple conditions. All in all, the current workaround of calling all boundary condition method might be redundant but at least works correctly. The only alternative I can think of would be to require all conditions to be enforced after the operation, remove the applyBefore... methods and be done with it.

Enough with boundary conditions for now. With the pieces we have so far, we can take a payoff at $$t=T$$ and evolve it back step by step to $$t=0$$—if nothing happens in between, that is. But that's not often the case, which is my cue for the next post.

## Monday, July 6, 2015

Hello everybody.

A very short post for a public service announcement. The folks at MoneyScience have organized a Google Hangout with me on next Monday, July 13th. Details are at this link.

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

## Monday, June 29, 2015

### A quick look at the QuantLib 1.6 release

Hello everybody.

Unlike in Westeros, summer is coming.  My posting schedule will be more relaxed, as usual, so don't hold your breath waiting for the next posts (hint: subscribe instead, and you'll be notified whenever I post).

A bit of news: the guys at Quants Hub are happy of the success of their new pricing structure, and they're extending it to July 31st. This means that, for the whole next month, you can still get my video workshop for £99.

However, the bigger piece of news is that last Tuesday I've released QuantLib 1.6; you can download it at this link if you haven't already. Thanks to all those who contributed. In this post, I'll share some information on the release (like I did last time).

I'm happy to report that we managed to decrease the distance between releases. This one comes out four and a half months after QuantLib 1.5, which is much better than the one-year hiatus between the 1.4 and 1.5 releases. However, this doesn't mean that the 1.6 release was a small one: it contains 65 issues and pull requests (follow the link to see their list and examine the changes brought by any one of them). A couple of quick git commands also shows a grand total of 382 commits by 16 contributors (git outputs 18 because it includes a couple of duplicates it doesn't know about).

As usual, the list doesn't include contributors that reported bugs, suggested changes, or even contributed patches but didn't actually author the relevant commits. Their names are in the list of changes for the 1.6 release; apologies if I missed anyone.

As I write, there's 27 open pull requests already in the queue for next release. Stay tuned.

Again, thanks to all those who contributed to the 1.6 release!

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

## Monday, June 22, 2015

### Chapter 8, part 2 of n: evolution schemes

Welcome back.

This week, the second part of the series on the finite-differences framework that started in the last post. And also this week, I'll be releasing QuantLib 1.6. Stay tuned.

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

### Evolution schemes

In a partial differential equation $${\partial f}/{\partial t} = L f$$, the operators in the previous subsection provide a discretization of the derivatives on the right-hand side. Evolution schemes discretize the time derivative on the left-hand side instead.

As you know, finite-difference methods in finance start from a known state $$\mathbf{f}(T)$$ at maturity (the payoff of the derivative) and evolve backwards to today's date. At each step, we need to evaluate $$\mathbf{f}(t)$$ based on $$\mathbf{f}(t+\Delta t)$$.

The simplest way is to approximate the equation above as
$$\frac{\mathbf{f}(t+\Delta t) - \mathbf{f}(t)}{\Delta t} = L \cdot \mathbf{f}(t+\Delta t)$$
which simplifies to $$\mathbf{f}(t) = \mathbf{f}(t+\Delta t) - \Delta t \cdot L \cdot \mathbf{f}(t+\Delta t)$$, or $$\mathbf{f}(t) = (I - \Delta t \cdot L) \cdot \mathbf{f}(t+\Delta t)$$. It's a simple matrix multiplication, and updating the function array translates into
    a = (I - dt*L).applyTo(a);

(in actual code, of course, the I + dt*L operator would be precalculated and stored across time steps).

A slightly less simple way is to write
$$\frac{\mathbf{f}(t+\Delta t) - \mathbf{f}(t)}{\Delta t} = L \cdot \mathbf{f}(t)$$
which results into $$(I + \Delta t \cdot L) \cdot \mathbf{f}(t) = \mathbf{f}(t + \Delta t)$$. This is a more complex step, since it requires solving a linear system (or inverting the operator, which is the same); but given the interface of our operators, the update translates equally simply into
    a = (I + dt*L).solveFor(a);

Now, this is where using tridiagonal operators pays off. Let me oversimplify, so I can keep a long story short: schemes of the first kind (called explicit schemes) are simpler but less stable. Implicit schemes, instead (those of the second kind), are more stable and can be used with larger time steps. If the operator were a generic one, the price for an implicit step would be $$O(N^3)$$ in the size of the grid, making it unpractical when compared to the $$O(N^2)$$ cost for an explicit step. For tridiagonal operator, though, both methods cost $$O(N)$$, which allows us to reap the benefits of explicit steps with no additional price.

Going back to the library: the two schemes above (called explicit Euler and implicit Euler scheme, respectively) are both available. They are two particular specializations of a generic MixedScheme class template, shown in the listing that follows, which models the discretization scheme
$$\frac{\mathbf{f}(t+\Delta t) - \mathbf{f}(t)}{\Delta t} = L \cdot \left[(1-\theta) \cdot \mathbf{f}(t+\Delta t) + \theta \cdot \mathbf{f}(t) \right]$$
where an implicit and explicit step are mixed. As you can see from the formula, $$\theta = 0$$ gives the explicit Euler scheme and $$\theta = 1$$ gives implicit Euler; any other value in the middle can be used (for instance, $$1/2$$ gives the Crank-Nicolson scheme).
    template <class Operator>
class MixedScheme  {
public:
typedef OperatorTraits<Operator> traits;
typedef typename traits::operator_type operator_type;
typedef typename traits::array_type array_type;
typedef typename traits::bc_set bc_set;
typedef typename traits::condition_type condition_type;

MixedScheme(const operator_type& L,
Real theta,
const bc_set& bcs)
: L_(L), I_(operator_type::identity(L.size())),
dt_(0.0), theta_(theta) , bcs_(bcs) {}

void setStep(Time dt) {
dt_ = dt;
if (theta_!=1.0) // there is an explicit part
explicitPart_ = I_-((1.0-theta_) * dt_)*L_;
if (theta_!=0.0) // there is an implicit part
implicitPart_ = I_+(theta_ * dt_)*L_;
}
void step(array_type& a, Time t) {
if (theta_!=1.0) {
a = explicitPart_.applyTo(a);
if (theta_!=0.0) {
implicitPart_.solveFor(a, a);
}
protected:
operator_type L_, I_, explicitPart_, implicitPart_;
Time dt_;
Real theta_;
bc_set bcs_;
};

template <class Operator>
class ExplicitEuler : public MixedScheme<Operator> {
public:
// typedefs, not shown
ExplicitEuler(const operator_type& L,
const bc_set> >& bcs)
: MixedScheme<Operator>(L, 0.0, bcs) {}
};

As I mentioned before, schemes take the type of their differential operators as a template argument. The MixedScheme class does that as well, and extracts a number of types from the operator class by means of the OperatorTraits template (that I'll describe later on).

The constructor of the MixedScheme class takes the differential operator $$L$$, the value of $$\theta$$, and a set of boundary conditions (more on that later) and stores them. Another method, setStep, stores the $$\Delta t$$ to be used as time step and precomputes the operators $$I - (1-\theta) \cdot \Delta t \cdot L$$ and $$I + \theta \cdot \Delta t \cdot L$$ on which the applyTo and solveFor methods will be called; the two special cases $$\theta = 0$$ and $$\theta = 1$$ are checked in order to avoid unnecessary computations. This is done in a separate method in case we needed to change the time step in the middle of a simulation; by being able to reset the time step, we can reuse the same scheme instance.

Finally, the step method performs the actual evolution of the function array. The implementation shown in the listing is just a basic sketch that I'll be improving in the next subsections; here, it just performs the explicit and implicit part of the scheme in this order (again, avoiding unnecessary operations).

The few schemes provided by the library (such as the ExplicitEuler, shown in the listing) are all specializations of the MixedScheme class (I'm just talking about the old framework here; the new framework has different ones). Again, there's no base class for an evolution scheme; other parts of the framework will take them as a template argument and expect them to provide the setStep and step methods.

## Monday, June 15, 2015

### Chapter 8, part 1 of n: the finite-differences framework.

Welcome back.

As I write, release candidates for QuantLib 1.6 have been out for a week or two (if you haven't heard about this already, you might want to subscribe to the QuantLib mailing list). If you have a few cycles to spare, please download it, give it a try and report any problems to the mailing list. The official 1.6 release should be out in a week or two. (Brace yourself for new developments afterwards. Just saying.)

This week, some initial content from chapter 8 of my book. Yes, I finally sat down to write it. The first couple subsections were pushed to the Leanpub edition a couple of weeks ago, with more coming during the summer (both there and here).

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

### The finite-differences framework

The last framework I'll be writing about was, in fact, the first to appear in the library. Of course, it was not a framework back then; it was just a bunch of code hard-coded inside a single option pricer.

Since then, the reusable parts of the code were reorganized and grew in complexity—probably to the point where it's hard to find one's way inside the architecture. A few years ago, a new version of the framework was contributed adding 2-dimensional finite-difference models (which were absent in the first version) and sporting a more modular architecture.

The two versions of the framework still coexist in the library. In this chapter, I'll describe both of them—while hoping for the new one to replace the first some day.

### The old framework

As I mentioned in a previous post, our master plan to create a common framework for trees and finite-differences models never came to fruition. Instead, the finite-differences code stayed closer to its underlying mathematics and implemented concepts such as differential operators, boundary conditions and so on. I'll tackle one of them in each of the following subsections.

#### Differential operators

If I were to start at the very beginning, I'd describe the Array class here; but they behave as you can imagine, so you'll forgive me if I leave its description to a later post and just note that it's used here to discretize the function $$f(x)$$ on a grid $$\{ x_0 \ldots x_{N-1} \}$$ as the array $$\mathbf{f}$$, with $$\mathbf{f}_i = f(x_i)$$ for any $$i$$.

Onwards to operators. Now, I'm probably preaching to the choir, so I'll keep the math to a minimum to avoid telling you lots of stuff you know already. (And if you don't, there's no way I can tell you all you need. You can refer to any of the several books on finite differences; for instance, Duffy, 2006 [1]).

In short: a differential operator transforms a function $$f(x)$$ in one of its derivatives, say, $$f^{\prime}(x)$$. The discretized version transforms $$\mathbf{f}$$ into $$\mathbf{f^{\prime}}$$, and since differentiation is linear, it can be written as a matrix. In general, the operator doesn't give the exact discretization of the derivative but just an approximation; that is, $$\mathbf{f^{\prime}}_i = f^{\prime}(x_i) + \epsilon_i$$ where the error terms can be reduced by decreasing the spacing of the grid.

QuantLib provides differential operators for the first derivative $$\partial/\partial x$$ and the second derivative $$\partial^2/\partial x^2$$. They are called $$D_0$$ and $$D_+D_-$$, respectively, and are defined as:
$$\frac{\partial f}{\partial x}(x_i) \approx D_0 \mathbf{f}_i = \frac{\mathbf{f}_{i+1}-\mathbf{f}_{i-1}}{2h} \qquad \frac{\partial^2 f}{\partial x^2}(x_i) \approx D_+D_- \mathbf{f}_i = \frac{\mathbf{f}_{i+1}-2\mathbf{f}_i+\mathbf{f}_{i-1}}{h^2}$$
where $$h = x_i - x_{i-1}$$ (we're assuming that the grid is evenly spaced). Taylor expansion shows that the error is $$o(h^2)$$ for both of them.

As you can see, the value of the derivative at any given index $$i$$ (except the first and last, that I'll ignore now) only depends from the values of the function at the same index, the one that precedes it, and the one that follows. This makes $$D_0$$ and friends tridiagonal: the structure of such operators can be sketched as
$$\left( \begin{array}{ccccccc} d_o & u_0 & & & & & \\ l_0 & d_1 & u_1 & & & & \\ & l_1 & d_2 & u_2 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & l_{n-4} & d_{n-3} & u_{n-3} & \\ & & & & l_{n-3} & d_{n-2} & u_{n-2} \\ & & & & & l_{n-2} & d_{n-1} \end{array} \right)$$

They have a number of desirable properties, including the fact that a linear combination of tridiagonal operators (such as the one found, say, in the Black-Scholes-Merton formula) is still tridiagonal; thus, instead of using a generic Matrix class (which does exist, but is not used here), we made them instances of a specific class called TridiagonalOperator and shown in the following listing.
    class TridiagonalOperator {
Size n_;
Array diagonal_, lowerDiagonal_, upperDiagonal_;

public:
typedef Array array_type;

explicit TridiagonalOperator(Size size = 0);
TridiagonalOperator(const Array& low,
const Array& mid,
const Array& high);

static Disposable<TridiagonalOperator>
identity(Size size);

Size size() const;
const Array& lowerDiagonal() const;
const Array& diagonal() const;
const Array& upperDiagonal() const;

void setFirstRow(Real, Real);
void setMidRow(Size, Real, Real, Real);
void setMidRows(Real, Real, Real);
void setLastRow(Real, Real);

Disposable<Array> applyTo(const Array& v) const;
Disposable<Array> solveFor(const Array& rhs) const;
void solveFor(const Array& rhs, Array& result) const;

private:
friend Disposable<TridiagonalOperator>
operator+(const TridiagonalOperator&,
const TridiagonalOperator&);
friend Disposable<TridiagonalOperator>
operator*(Real, const TridiagonalOperator&);
// other operators
};

As you can see from the private members, the representation of the operator reflects its tridiagonal structure: instead of storing all the mostly null elements of the matrix, we just store three arrays corresponding to the diagonal, lower diagonal, and upper diagonal; that is, the $$d_i$$, $$l_i$$ and $$u_i$$ of the earlier sketch.

The constructors and the first bunch of methods deal with building and inspecting the operator. The first constructor sets the size of the matrix (that is, of the three underlying arrays) and trusts the values to be filled later, by means of the provided setters: referring to the figure again, the setFirstRow method sets $$d_0$$ and $$u_0$$, the setLastRow method sets $$l_{n-2}$$ and $$d_{n-1}$$, setMidRow sets $$l_{i-1}$$, $$d_i$$ and $$u_i$$ for a given $$i$$, and the convenience method setMidRows sets them for all $$i$$. The second constructor sets all three arrays in one fell swoop; and a factory method, identity, works as a third constructor by building and returning an operator with $$1$$ on the diagonal and $$0$$ elsewhere. The inspectors do the obvious thing and return each the corresponding data member. Oftentimes, the returned values are wrapped in a Disposable class template, in a possibly misguided attempt to avoid a few copies. Again, details will be in a future post.

The applyTo and solveFor methods are the meat of the interface. Given a tridiagonal operator L and an array u, L.applyTo(u) returns the result $$L u$$ of applying $$L$$ to $$u$$, while L.solveFor(u) returns the array $$v$$ such that $$u = L v$$ (The overloaded call L.solveFor(u,v) does the same while avoiding an allocation). Both operations will be used later, and for a tridiagonal operator they both have the nice property to be $$O(N)$$ in time, that is, to only take time proportional to the size $$N$$ of the operator. The algorithms were taken from the classic Numerical Recipes in C [2]. Not the code, of course: besides not being free, the original code had a number of endearing old conventions (such as numbering arrays from 1, like in Downton Abbey) that wouldn't have matched the rest of the library.

Finally, the class defines a few operators. They make it possible to write numerical expressions such as 2*L1 + L2 that might be used to compose tridiagonal operators. To this effect, the library also provides a few building blocks such as the $$D_0$$ and $$D_+D_-$$ operators I mentioned earlier; for instance, given the Black-Scholes equation, written in operator form as
$$\frac{\partial f}{\partial t} = \left[ - \left(r - q - \frac{\sigma^2}{2} \right) \frac{\partial}{\partial x} - \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2} + r I \right] f \equiv L_{BS} f$$
one could define the corresponding $$L_{BS}$$ operator as
    TridiagonalOperator L_BS = -(r-q-sigma*sigma/2)*DZero(N)
-(sigma*sigma/2)*DPlusDMinus(N)
-r*TridiagonalOperator::identity(N);

In this case, though, we didn't eat our own dog food. The library does define a Black-Scholes operator, but not through composition as above.

A final note: while they're probably the most performant, tridiagonal operators are not the only ones that could be used in a finite-difference scheme. However, there's no base class for such operators. The interface that they should implement (that is, applyTo and solveFor is only defined implicitly, through its use in the template classes that we'll see in the next subsection. In hindsight, defining such a base class wouldn't have degraded performance; the cost of a virtual call to applyTo or solveFor would have been negligible, compared to the cost of the actual calculation. But we were eagerly reading about template numerical techniques at that time; and like the two guys in the Wondermark strip, we decided to jingle all the way.

#### Bibliography

[1] D. J. Duffy, Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach. John Wiley and Sons, 2006.
[2] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd edition. Cambridge University Press, 1992.