## 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.