`PricingEngine`

and the `GenericPricingEngine`

class. These engines provide an interface and some basic functions but the actual calculations have not been implemented. One of the simplest pricing engines that actually performs some calculations is the `DiscountingBondEngine`

. This engine uses a discount curve to discount all the cash flows. In fact, the implementation of simply calls the `CashFlows::NPV`

method to evaluate the NPV and the settlement value.
The public interface of the `DiscountingBondEngine`

looks as follows.

class DiscountingBondEngine : public Bond::engine { public: DiscountingBondEngine( const Handle<YieldTermStructure>& discountCurve = Handle<YieldTermStructure>(), boost::optional<bool> includeSettlementDateFlows = boost::none); void calculate() const; Handle<YieldTermStructure> discountCurve(); };

The constructor takes a handle to a `YieldTermStructure`

which contains the discount curve. An optional boolean can be specified to indicate whether cash flows on the settlement date should be included in the calculations. If the flag is not specified, the global value of the flag from the Settings singleton will be used. One more thing to be noted is that the `DiscountingBondEngine`

inherits from `Bond::engine`

and not directly from `PricingEngine`

. The reason is that the `Bond`

class checks the type of the pricing engine that is passed to it. If the eninge cannot be cast to a `Bond::engine`

an exception is thrown. The definition of the `Bond::engine`

class is simple.

class Bond::engine : public GenericEngine<Bond::arguments, Bond::results> {};

This engine does not add any functionality but simply defines the arguments and results that the pricing engine should use. As explained above, any pricing engine that is to be used with the Bond instrument class MUST inherit from `Bond::engine`

. The reason for this is that the Bond class expects the type of arguments to be of type `Bond::arguments`

and the results to be of type `Bond::results`

. More about arguments and result types below.

Remember in the discussion of the Instrument class that arguments should be set inside the virtual `setupArguments`

method. The parameter to the method is a pointer to `PricingEngine::arguments`

. The pricing engine should supply the actual arguments object via the `getArguments`

method. By virtue of the `GenericEngine`

template class we can ensure that the arguments will, in fact, always be of type `Bond::arguments`

.

class Bond::arguments : public PricingEngine::arguments { public: Date settlementDate; Leg cashflows; Calendar calendar; void validate() const; };

The arguments include a settlement date, a leg of cash flows and a calendar that is used to adjust the cash flow dates. The validate method ensures that a settlement date is specified and that there is at least one cash flow in the leg. When an object of this type is passed to the `setupArguments`

method of the `Bond`

class, it is filled with the values taken from the bond.

The `Instrument`

class defines, and already implements the `fetchResults`

method. As I stated __here__ the `Instrument`

class fetches the NPV, the error estimate and any additional results from the pricing engine. It also sets the `valuationDate`

of the instrument. In order to do this, the results type returned by the `PricingEngine`

has to be of the nested type `Instrument::results`

which extends `PricingEngine::results`

.

class Instrument::results : public virtual PricingEngine::results { public: void reset(); Real value; Real errorEstimate; Date valuationDate; std::map<std::string,boost::any> additionalResults; };

This class is filled by the pricing engine and delivers the results to the `Instrument`

class which then reads from it in the `fetchResults`

method. In addition to the fields defined in the `Instrument::results`

class, the `Bond`

class requires one additional value to be calculated by the pricing engine, namely the settlement value. For this reason `Bond`

defines its own results type. Because `Bond`

extends the `Instrument`

class the result type defined by `Bond`

also has to extend `Instrument::results`

.

class Bond::results : public Instrument::results { public: void reset(); Real settlementValue; };

As stated above, the `Bond::engine`

ensures that it can only take argument types and return result types that extend from `Bond::arguments`

and `Bond::results`

by extending `GenericEngine<Bond::arguments, Bond::results>`

.

All of what I described in the last two sections happens behind the scenes and is only interesting if you intend to write your own pricing engines. If you want to simply use the existing pricing engines that QuantLib provides you simply need to pass a pricing engine to the `setPricingEngine`

method of the instrument. The only thing that you need to ensure is that you are using a pricing engine class that is compatible with your instrument. Note again that this compatibility is not checked during compile time. Instead an exception is thrown during runtime if a non-compatible pricing engine is used.

If you are interested in the internals of QuantLib and you are thinking about implementing your own pricing engines or even you own instruments here is a little diagram that might help you. I have extracted the UML diagram of some of the classes that are related to the `DiscountingBondEngine`

. In this hierarchy the `PricingEngine`

and its nested classes `arguments`

and `resuls`

stand at the top. In this diagram I have omitted the properties and methods of the classes involved for clarity.

QuantLib provides a number of functions that calculate the NPV of a series of cash flows. These are defined as static methods inside the `CashFlows`

class. The simplest calculation of the NPV will take a constant interest rate over the period of the cash flows.

static Real npv(const Leg& leg, const InterestRate& yield, bool includeSettlementDateFlows, Date settlementDate = Date(), Date npvDate = Date()); static Real npv(const Leg& leg, Rate yield, const DayCounter& dayCounter, Compounding compounding, Frequency frequency, bool includeSettlementDateFlows, Date settlementDate = Date(), Date npvDate = Date());

The first function takes a `Leg`

of cash flows and an interest rate that defines the yield. The interest rate is passed using the InterestRate class. In short, `InterestRate `

contains information about interest rates such as the rate itself, the compounding rule, the frequency and a day counting convention. Note that this allows the `InterestRate`

class to calculate a discount factor for a given time interval. The settlement date can be passed as fourth argument. Only cash flows that occur on or after the settlement date will contribute to the NPV. The inclusion of cash flows on the settlement date is governed by the flag `includeSettlementDateFlows`

. As always the parameters are simply passed to the `hasOccured`

method of the cash flow to determine if the cash flow has already occurred. If the settlement date is not set the global evaluation date will be used. The last parameter can be used to specify a different date to which the cash flows are discount back to. This means that the discounting factor is calculated using `npvDate`

. If `npvDate`

is not set, the settlement date will be taken.

Note that the interest rate is calculated by compounding the factor calculated for the cash flow intervals. What this means is that the discount rate is calculated using the compounding rules defined in the `InteresRate yield`

object for the intervals between two cash flows. Let’s say that \(d_i\) and \(d_{i-1}\) are the dates of two successive cash flows, then the interval discount factor is calculated as a function of these two dates.

$$\mathrm{Discount Factor}_i = DInt_{i, i-1}(d_i, d_{i-1})$$

Here \(DInt\) is the discount rate calculated by the `InterestRate`

object and might use simple, continuous compounding or use a compounding discount factor with a frequency that does not necessarily match the frequency of the cash flows. To get the discount factor of a cash flow these interval discount factors are compounded by multiplying them together.

$$D_i = \prod_{i=1}^{N} DInt_{i, i-1}(d_i, d_{i-1})$$

The final net present value is then calculated using the individual discount factors.

$$\mathrm{NPV} = \sum_{i=1}^{N} D_i \times R_i$$

In many realistic cases, one won’t have a constant interest rate over the period of the cash flows. In this case the NPV of a `Leg`

has to be calculated using a yield curve. In QuantLib yield curves are represented by `YieldTermStructure`

.The `CashFlows`

class provides a static method for evaluating the NPV in this case.

Real npv(const Leg& leg, const YieldTermStructure& discountCurve, bool includeSettlementDateFlows, Date settlementDate = Date(), Date npvDate = Date());

As with the two methods described above, `includeSettlementDateFlows`

determines if cash flows on the settlement date are taken into account for the calculation. The date `settlementDate`

is the settlement date which enters into the `hasOccured`

method of the cash flows. Again, this means only cash flows on or after the settlement date are included in the calculation. The `npvDate`

is a date that the NPV is discounted back to. To understand what the function calculates, assume that \(D(d_i)\) is the discount factor returned by the yield curve for a given date \(d_i\). This discount factor is simply the return value of the `discount`

method of `YieldTermStructure`

. The NPV is then calculated using the formula

$$\mathrm{NPV} = \frac{\sum_{i=1}^{N} D(d_i) \times R_i}{D(\mathrm{npvDate})}$$

The last function for calculating the NPV takes the zero spread of the cash flows into account.

Real npv(const Leg& leg, const boost::shared_ptr<YieldTermStructure>& discount, Spread zSpread, const DayCounter& dayCounter, Compounding compounding, Frequency frequency, bool includeSettlementDateFlows, Date settlementDate = Date(), Date npvDate = Date());

This function uses the zero spread, given by `zSpread`

, the day counter, the compounding rule and the frequency to construct a modified yield curve. This modified yield curve is simply the original curve passed to the function with an added zero spread. The modified curve is then passed to the previous function for calculating the NPV.

`CashFlow`

class, typedef std::vector< boost::shared_ptr<CashFlow> > Leg;

The functions for analysing cash flows are, somewhat confusingly, defined as static member functions of a class called `CashFlows`

. Notice the plural **s** at the end and don’t confuse it with the `CashFlow`

class that holds information about a single cash flow event. Many methods assume that the cash flows in a leg are sorted by their date as given by `Event::date()`

. This will not be pointed out individually for the methods described below.

In this post I will be talking about some basic methods defined inside the `CashFlows`

. These methods can be used to obtain dates or amounts of specific cash flows. In subsequent posts will focus on the methods that evaluate the cash flows and calculate data such as the NPV.

You can find the start date and end date of a `Leg`

using the `startDate`

and `endDate`

methods. The method `isExpired`

will return true when the leg is expired and false if it isn’t.

static Date startDate(const Leg& leg); static Date maturityDate(const Leg& leg); static bool isExpired(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date());

`startDate`

will find and return the earliest date in `leg`

. For most cash flows, the dates are taken from the `date`

method of the `Event`

class. For `Coupon`

cash flows, on the other hand, the start date of the accrual period is obtained by calling the `accrualStartDate`

method. Similarly `maturityDate`

will return the latest cash flow in the leg. For `Coupon`

cash flows the `accrualEndDate`

is taken to be the relevant date, for all other cash flows the normal date returned by `date`

is used.

The method `isExpired`

will query all cash flows if they have occurred. If all cash flows return `true`

then the leg is expired and `isExpired`

returns true. If one of the cash flows returns `false`

when queried via the `hasOccured`

method then `isExpired`

will return `false`

. The `settlementDate`

and `includeSettlementDateFlows`

parameters are simply passed on to the `hasOccured`

method in the `Event`

class. If the settlement date is not specified then the `Settings::evaluationDate`

will be used. An empty leg is considered expired.

static Leg::const_reverse_iterator previousCashFlow(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static Leg::const_iterator nextCashFlow(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static Date previousCashFlowDate(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static Date nextCashFlowDate(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static Real previousCashFlowAmount(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static Real nextCashFlowAmount(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date());

The method `previousCashFlow`

will return the last cash flow on or before a given date while `nextCashFlow`

will return the next cash flow after a given date. The methods use the `hasOccured`

method of the `Event`

class to decide if a cash flow has occurred. As with `isExpired`

, the `settlementDate`

and `includeSettlementDateFlows`

parameters are passed on to `hasOccured`

and an empty settlement date will mean that the global evaluation date is used.

`previousCashFlow`

will return a reverse iterator pointing to the last cash flow that has occurred while `nextCashFlow`

will return a forward iterator to the first cash flow that has not occurred. If no cash flow was found that matches the critetria then `previousCashFlow`

will return `leg.rend()`

while `nextCashFlow`

will return `leg.end()`

. The advantage of returning a reverse iterator from `previousCashFlow`

is that you can easily iterate over all cash flows that have already occurred.

The other four methods are convenience methods that return the date or amount of the cash previous or next cash flow. `previousCashFlowDate`

calls `previousCashFlow`

and returns the date of the cash flow, and similarly `nextCashFlowDate`

calls `nextCashFlow`

and returns the date. In both cases a null `Date`

is returned if no cash flow was found.

`previousCashFlowAmount`

and `nextCashFlowAmount`

will return the amount of the previous or next cash flow. Again the cash flows are found by calling `previousCashFlow`

or `nextCashFlow`

. If there are multiple cash flows that occur at the same date then all their amounts are added together. If no cash flow is found a value of zero is returned.

QuantLib defines many different types of cash flows and, using polymorphism, a `Leg`

array can hold objects of different cash flow classes. However, in some cases one might be interested only in cash flows of a single type. Specifically, there are many cases in which only the `Coupon`

cash flows should be taken into account. The `CashFlows`

class defines a number of methods that only consider coupons within a cash flow and return `Coupon`

related information. In this section I will discuss only the functions that return basic information about coupon dates and periods. The functions that perform calculations on coupons are dealt with in a future post.

static Date accrualStartDate(const Leg& leg, bool includeSettlementDateFlows, Date settlDate = Date()); static Date accrualEndDate(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static Date referencePeriodStart(const Leg& leg, bool includeSettlementDateFlows, Date settlDate = Date()); static Date referencePeriodEnd(const Leg& leg, bool includeSettlementDateFlows, Date settlDate = Date()); static Time accrualPeriod(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static BigInteger accrualDays(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static Time accruedPeriod(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static BigInteger accruedDays(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static Real accruedAmount(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date());

If you look at the documentation of coupons, you will notice that these methods are direct mirror images of the methods defined by the `Coupon`

class. All these methods will find the first `Coupon`

cash flow in the leg that has not yet occurred. Then they call the respective method on that coupon.

For example the method `accruedPeriod`

defined within `CashFlows`

will first find the first `Coupon`

object in leg that occurred after `settlementDate`

. Once a coupon, let’s call it `cp`

, is found the method will return `cp.accruedPeriod`

. If no coupon is found a sensible default value will be returned, in this case a period of 0.0. In all methods that return a numerical value the default is zero, for methods that return a date the default is an empty `Date`

object. As in `nextCashFlow`

, if the settlement date is not specified then the `Settings::evaluationDate`

will be used. The `hasOccurred`

method is called, with the settlement date and the `includeSettlementDate`

parameter, to find the cash flow that has not yet occurred.

Two Coupon related function behave slightly differently to the functions above.

static Rate previousCouponRate(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date()); static Rate nextCouponRate(const Leg& leg, bool includeSettlementDateFlows, Date settlementDate = Date());

The function will return the coupon rate of the previous or the next coupon in the leg, with respect to the settlement date. In some cases the leg might contain multiple coupons for the same date. If this is the case then an aggregate rate will be calculated for all coupons that occur on the same date. This aggregate rate is simply the sum of all the rates of the coupons on that date. However, for the calculation to be successful all coupons on that date must have the same nominal value, accrual period and must be using the same day counter.

]]>`boost/random`

package, you will find a wealth of random number generators with different properties. In some cases you might want to test if your results depend on the specific implementation of the random number generator and whether some other generator might give better, or more dependable results. In this post I will write two wrapper classes that will allow any Boost random number generator to be used with QuantLib.
To understand the differences between the Boost generators and the generators of QuantLib, I will first analyse the anatomy of the Boost random number generators. Boost knows two types of random number generators, integral and floating point generators. As the name suggests, integral generators return integral values and floating point generators return floating point values. Let’s look at the integral generators first. The integral type that a generator returns can have different bit length. For example the `mt19937_64`

is a 64 bit version of the Mersenne twister with periodicity of \(2^{19937}-1\). Every generator class in the Boost library defines the `result_type`

as the type that the generator returns. For example `rand48`

contains the line

typedef boost::uint32_t result_type;

All Boost random number generators override the function operator to access the next random number.

result_type operator()();

This will generate and return the next number in the sequence. Random number generators all have a default constructor and a constructor that takes a seed. For all integral engines the seed type is the same as the result type. In addition there is a `seed`

method that re-initialises a random number generator with a seed. For example `rand48`

has the following constructors.

rand48(); rand48(result_type seed);

And the `seed`

method is defined as follows.

seed(result_type seed);

There is also a `seed`

function without arguments but we will only discuss the methods that are needed to create our generic wrapper class.

The aim is to write a generic wrapper around all of Boost’s random number generators. The idea is to follow the pattern used in the `Ranlux3UniformRng`

wrapper for Boost’s `ranlux64_3_01`

RNG. The difference will be that this wrapper will be templated with the type of the Boost RNG. I will start by creating the `GenericBoostIntRNG`

class with the template argument RNG

template<class RNG> class GenericBoostIntRNG {

Here `RNG`

is the tye of the Boost random number generator. `GenericBoostIntRNG`

needs to keep an instance of `RNG`

as private member. Here we encounter the first little hurdle. Getting the next random number from a Boost RNG will modify the object. The reasoning is that the internal state of the RNG is changed and this manifests itself in the fact that the function operator is non-const. QuantLib, on the other hand, expects that getting the next random number does not modify the generator. The `next`

method is a `const`

method in all the QuantLib RNGs. The solution is to make the Boost RNG a mutable member of the wrapper.

private: mutable RNG generator;

Next we have to declare two local types. QuantLib demands the `sample_type`

to define the return type of the next method. On the other hand, we also need to use `result_type`

of `RNG`

in many places and should therefore declare a shortcut.

public: typedef Sample<Real> sample_type; typedef typename RNG::result_type result_type;

To make sure that we use `GenericBoostIntRNG`

only with integral random number generators we should make sure that `result_type`

is an integer type. This can be done with `BOOST_STATIC_ASSERT`

and the standard `numeric_limits`

type.

BOOST_STATIC_ASSERT (std::numeric_limits<result_type>::is_integer);

At this point we should be defining the constructor. This is not straightforward and will require some extra code because we might need to create an automatic seed. We will postpone this and first focus on the easier methods that generate the random numbers.

Real nextReal() const { return (Real(generator()) + 0.5) / ceil_value; }

This method creates a random integer number and convert it into a `Real`

value. The `ceil_value`

will be one greater than the maximum number returned by the random number generator. We define it as a static private member of the class.

private: static Real ceil_value;

The value is set by using the `max`

function of the Boost random number generator.

template<class RNG> Real GenericBoostIntRNG<RNG>::ceil_value = Real(RNG::max()) + 1.0;

With this definition, the call to generator in the `nextReal`

method will always return a value from `0`

to `ceil_value-1`

. This means that the value returned by `nextReal`

will always fall within the closed interval.

$$\left[ \frac{0.5}{ceil}, \frac{ceil-0.5}{ceil}\right]$$

Mathematically this interval will never include 0 or 1 but we have to be careful when stating that the numerical value will lie in the open interval \((0,1)\). Rounding errors might cause the upper value to round to 1. This will be the case when the `Real`

type has less bits that the integral type of the Boost RNG. If we stick to using 32 bit Boost RNGs and define Real as 64 bit `double`

than all is fine and the numbers will always be less than 1. In any case the number will never be larger than 1 and an occasional value of 1 should not cause a problem. QuantLib random number generators are expected to define a `next`

method that returns a `sample_type`

. This is quickly done.

public: sample_type next() const { return sample_type(nextReal(), 1.0); }

Now we come back to creating the constructor for our generic wrapper class.

public: GenericBoostIntRNG(seed_type seed = 0) : generator(getSeed(seed)) {}

This does not look complicated. We are initializing the Boost object generator with a seed that we obtain from a method called `getSeed`

. The `getSeed`

method will have to be implemented. If we specify a seed that is non-zero that we will simply pass it on. Otherwise we will use the QuantLib `SeedGenerator`

class to create a suitable seed for the generator object.

private: result_type getSeed(result_type seed) { unsigned long s = (seed != 0 ? seed : SeedGenerator::instance().get()); return (result_type) (s % RNG::max()); }

In the first line inside the method we assign the seed to an `unsigned long`

. The `SeedGenerator`

class in QuantLib is defined as a singleton. A call to `SeedGenerator::instance().get()`

will return an `unsigned long`

seed value based on the current time. Because `unsigned long`

might be longer than the `result_type`

of the Boost RNG we have to restrict the values by taking the modulus with the maximum value of the RNG.

With this the generic wrapper is complete. I have tested it by plugging various Boost random number generators into the EquityOption example of QuantLib. When using the Boost’s mt19937 generator, I get identical results to the QuantLib MersenneTwisterUniformRng.

Boost also specifies a number of floating point random number generators. For completeness I have also created a wrapper for using these within QuantLib. Using the experience gained from GenericBoostIntRNG the implementation of GenericBoostFloatRNG is relatively straightforward.

template<class RNG> class GenericBoostFloatRNG { private: mutable RNG generator; static Real ceil_value; public: typedef Sample<Real> sample_type; typedef typename RNG::result_type result_type; typedef boost::uint32_t seed_type; BOOST_STATIC_ASSERT (std::numeric_limits<result_type>::is_float); GenericBoostFloatRNG(seed_type seed = 0) : generator(getSeed(seed)) {} sample_type next() const { return sample_type(nextReal(), 1.0); } result_type nextFloat() const { return generator(); } private: seed_type getSeed(seed_type seed) { unsigned long s = (seed != 0 ? seed : SeedGenerator::instance().get()); return (seed_type) (s % 0xffffffff); } };

There are three differences between this class and the `GenericBoostIntRNG`

class. First, the `nextFloat`

method can return the result of the generator directly. Second, we need to define the `seed_type`

separately from the `result_type`

. The `seed_type`

will be integral while the `result_type`

will be a floating point type. Third, the restriction of the seed to values accepted by the Boost RNG looks slightly different. All floating point random number generators accept `uint32`

as seed type. This means we can simply restrict the seed by applying a bit mask.

I hope someone will find this useful. The code for these wrapper classes is freely available under the MIT licence and available on Github.

Follow the author

]]>`std::vector`

with numbers. While random number generators can easily be used to create a random sequence, some algorithms, notably the low discrepancy sequence generators, can only be used to generate sequences and not individual random numbers.
QuantLib defines four pseudo-random number generators with uniform distribution. These are the Mersenne Twister 19937, the L’Ecuyer generator with added Bays-Durham shuffle, Knuth’s random number generator and Luescher’s Luxury random number generator. The last generator in this list is not implemented inside QuantLib but is a wrapper for the random number generator implemented in the **Boost Random** library. There are two versions of the luxury random number generator, Ranlux 3 and Ranlux 4. The second generator is slower than the first but provides more “luxury” (see M. Luescher, Computer Physics Communications, **79** (1994) pp 100-110). The number generators are defined in the following classes.

class MersenneTwisterUniformRng; class LecuyerUniformRng; class KnuthUniformRng; class Ranlux3UniformRng; class Ranlux4UniformRng;

All of these generators have constructors that take an optional initial seed as argument.

explicit MersenneTwisterUniformRng(unsigned long seed = 0); explicit LecuyerUniformRng(long seed = 0); explicit KnuthUniformRng(long seed = 0); explicit Ranlux3UniformRng(Size seed = 19780503u); explicit Ranlux4UniformRng(Size seed = 19780503u);

As you can see, different generators use slightly different integral types as seeds. The types depend on the internals of the generators. The default seed for the Mersenne Twister, L’Ecuyer and Knuth generators will cause a seed to be generated based on the system clock. The seed is calculated from the return value of the `time()`

function using the Mersenne Twister. All this is carried out inside the `SeedGenerator`

class.

While the details of the seeding might not be important to you there is one thing to take away from this. The seed is calculated in a deterministic way from the result of the `time()`

function but this function has a granularity of one second. If you run the same code multiple times during one second, the random number generator is likely to produce the same values.

`MersenneTwisterUniformRng`

provides another constructor that takes a vector of seeds which will initialise the internal states of the Mersenne Twister.

explicit MersenneTwisterUniformRng( const std::vector& seeds);

Each of the four random number generators has a method to get the next random number. The random number returned is a `Real`

number between, but not including, 0 and 1.

Sample<Real> next() const;

The return type is actually not just a simple `Real`

but a `Real`

value together with a weight wrapped inside the `Sample`

struct.

template <class T> struct Sample { typedef T value_type; Sample(const T& value, Real weight) : value(value), weight(weight) {} T value; Real weight; };

This struct is used throughout the QuantLib random number module to represent random numbers. For all samples returned by the random number generators, the `weight`

property is set to 1.0. To illustrate how to use the random number generators to create and print a sequence of random numbers, let’s look at a simple example.

MersenneTwisterUniformRng generator(12345); for (int i=0; i<100; ++i) { std::cout << generator.next().value << std::endl; }

The code above will create a Mersenne Twister and seed it with the number 12345. It will then create and print out 100 numbers.

You can easily create a random sequence from a random number generator with the generic `RandomSequenceGenerator`

. This class is templated with a random number generator class so it can use any generator to create the sequence.

template<class RNG> class RandomSequenceGenerator;

The class defines a type called `sample_type`

that defines the return type of the sequence generator. In general all sequence generators define the sample_type in the following way.

typedef Sample<std::vector<Real> > sample_type;

The generic random sequence generator has two constructors. The first argument of both constructors is the dimension, i.e. the length of the sequence that the generator will return. The second parameter is either a random number generator of the type specified by `RNG`

or a seed number that is used to initialise the random number generator. The seed is optional and defaults to zero if not specified.

RandomSequenceGenerator::RandomSequenceGenerator (Size dimensionality, const RNG& rng); RandomSequenceGenerator::RandomSequenceGenerator (Size dimensionality, BigNatural seed = 0);

Two methods allow access to the sequence.

const sample_type& nextSequence() const; const sample_type& lastSequence() const;

The method `nextSequence`

generates a new sequence with a length specified by dimensionality and returns a reference to it. Repeated access to the last sequence, without generating it again, is granted by `lastSequence`

. The following example shows how the sequence generator might be used.

RandomSequenceGenerator<MersenneTwisterUniformRng> generator(100, 12345); std::vector<Real> values = generator.nextSequence().value; for (int i=0; i<values.size(); ++i) { std::cout << values[i] << std::endl; }

All the previous random number and random sequence generators produce uniformly distributed random numbers. In many cases random numbers with different probability distributions are needed. The most common distribution is the Gaussian distributions. QuantLib provides two classes that generate Gaussian distributed random numbers from uniformly distributed random numbers. The `BoxMullerGaussianRng`

uses the Box Muller transform that generates pairs of normally distributed values from pairs of uniformly distributed values. Because the Box Muller transform needs to evaluate the square root and the logarithm for every number produced it can be slow for some applications. An alternative is the application of the central limit theorem, implemented in `CLGaussianRng`

. This generator uses the fact that the sum of 12 uniformly random numbers is already a good approximation for a normal distributed random variable. Both classes are templated with the type of the underlying random number generator.

template <class RNG> class BoxMullerGaussianRng; template <class RNG> class CLGaussianRng;

Both classes are simply constructed by passing the random number generator to the constructor.

explicit BoxMullerGaussianRng::BoxMullerGaussianRng (const RNG& uniformGenerator); explicit CLGaussianRng::CLGaussianRng (const RNG& uniformGenerator);

In addition, both classes follow the same pattern as the uniform random number generators for accessing the numbers and define the `next`

method.

Sample<Real> next() const;

In this way, both Gaussian generators can be plugged into the `RandomSequenceGenerator`

in the same way as any of the uniform generators. The following example will illustrate this.

typedef BoxMullerGaussianRng<MersenneTwisterUniformRng> BoxMuller; RandomSequenceGenerator<BoxMuller> generator(BoxMuller(100, 12345)); std::vector<Real> values = generator.nextSequence().value; for (int i=0; i<values.size(); ++i) { std::cout << values[i] << std::endl; }

The code above creates a vector filled with 100 normally distributed values. Inside the loop the values are then printed out.

For distributions other than Gaussian distributions QuantLib defines a general class called `InverseCumulativeRng`

. The class creates random numbers with an arbitrary distribution from its inverse cumulative distribution. Assume that the desired probability distribution of the random numbers is given by the function \(f(x)\). Then the cumulative distribution function is then given by

$$F(x) = \int_{-\infty}^x f(x’) dx’$$

All that is needed is a function object that evaluates the inverse of the cumulative distribution \(F^{-1}(x)\). The declaration of the `InverseCumulativeRng`

class looks like this.

template <class RNG, class IC> class InverseCumulativeRng { public: typedef Sample<Real> sample_type; typedef RNG urng_type; explicit InverseCumulativeRng(const RNG& uniformGenerator); sample_type next() const; };

The class `IC`

must have a default constructor and must overload the function operator to return the inverse cumulative distribution.

Real IC::operator()(Real x);

The following example may be helpful in illustrating how to use the `InverseCumulativeRng`

class. Let’s assume that you want to create a random variable with a probability distribution given by the logistic distribution with mean zero and scale parameter 1,

$$f(x) = \frac{1}{\left(e^{x/2} + e^{-x/2}\right)}$$

Then the cumulative distribution function is given by

$$F(x) = \int_{-\infty}^x \frac{1}{\left(e^{x’/2} + e^{-x’/2}\right)} dx’ = \frac{1}{1+e^{-x}}$$

The inverse of the function can be found by solving \(y=F(x)\) for \(x\) which results in

$$F^{-1}(y) = \mathrm{ln}\left(\frac{y}{1-y}\right)$$

All we need to do is to create a class that evaluates this function.

class ICLogistics { public: Real operator()(Real y) { return log(y/(1-y)); } };

Now we can define the following random number generator.

typedef InverseCumulativeRNG< MersenneTwisterUniformRng, ICLogistics > LogisticsRNG;

In the example above we can now simply replace `BoxMuller`

with `LogisticsRNG`

to produce a sequence of random numbers with the logistics distribution.

Follow the author

]]>`TermStructure`

that defines methods common to any kind of term structure. It was pointed out that many classes inherit from the `TermStructure`

class, including a class called `YieldTermStructure`

. This class is, again, an abstract class and it defines the base for any yield curves. A yield curve is defined by a list of quotes together with a list of jump dates. These are passed to `YieldTermStructure`

in the constructor and are then stored internally. There are three different constructors which correspond to the three constructors defined in `TermStructure`

.
YieldTermStructure( const DayCounter& dc = DayCounter(), const std::vector<Handle<Quote> >& jumps = std::vector<Handle<Quote> >(), const std::vector<Date>& jumpDates = std::vector<Date>()); YieldTermStructure( const Date& referenceDate, const Calendar& cal = Calendar(), const DayCounter& dc = DayCounter(), const std::vector<Handle<Quote> >& jumps = std::vector<Handle<Quote> >(), const std::vector<Date>& jumpDates = std::vector<Date>()); YieldTermStructure(Natural settlementDays, const Calendar& cal, const DayCounter& dc = DayCounter(), const std::vector<Handle<Quote> >& jumps = std::vector<Handle<Quote> >(), const std::vector<Date>& jumpDates = std::vector<Date>());

All three constructors take two additional arguments, a `std::vector`

of `Quote`

handles and a `std::vector`

of jump dates. These two vectors define the data for the interest rate term structure. If the `jumpDates`

array is left empty, the jumps are assumed to occur annually at 31 Dec every year, starting in the year of the reference date. For the meaning of the remaining arguments and the difference between the constructors, see the documentation of the TermStructure class.

Two inspector methods give access to the jump dates.

const std::vector<Date>& jumpDates() const; const std::vector<Time>& jumpTimes() const;

The method `jumpDates`

returns vector of Date objects while `jumpTimes`

returns a vector of times. Recall that `Time`

is just a `Real`

number. Times are expressed in fractions of a year from the reference date and are evaluated using the day counter. They are calculated using the `timeFromReference`

method of `TermStructure`

.

`YieldTermStructure`

also overrides the `update`

method that it inherited from the Observer interface. Whenever the reference date changes, the jump times have to be recalculated.

Two methods return the discount factor from an arbitrary date to the reference date.

DiscountFactor discount(const Date& d, bool extrapolate = false) const; DiscountFactor discount(Time t, bool extrapolate = false) const;

The two methods differ in the way that the date is passed to them. The first version takes a date, while the second version takes a time measured in fractions of a year. This time should be calculated using the same day counting convention that is used in the `YieldTermStructure`

. In fact, the first method is just a convenience method which converts the `Date`

into `Time`

using `TermStrucure::timeFromReference`

and then calls the second version of the `discount`

method.

The discount is calculated by multiplying all the jumps that fall in the period from the reference date to the date passed to the discount() method and then multiplying this with the result of `discountImpl`

, discussed later. The `extrapolate`

flag indicates whether the discount factor should be extrapolated if it falls outside of the date range of the term structure.

The `discount`

method is used internally in a number of convenience methods. `zeroRate`

calculates the implied zero-coupon interest rate for a given time interval.

InterestRate zeroRate(const Date& d, const DayCounter& resultDayCounter, Compounding comp, Frequency freq = Annual, bool extrapolate = false) const; InterestRate zeroRate(Time t, Compounding comp, Frequency freq = Annual, bool extrapolate = false) const;

As before, the two methods differ in the way that the date is passed to them. The first version takes a specific date. This version allows passing a different day counter than the one used in the term structure. The `InterestRate`

that is returned uses the new day counter, compounding rule and frequency. The second version takes a time interval which has to be calculated using the same day counter as the one used in `YieldTermStructure`

. Internally, both methods first calculate the compound rate as the reciprocal of the discount rate, *a=1/d*. Then an `InterestRate`

is generated by calling `InterestRate::impliedRate()`

with the correct arguments. The extrapolate flag indicates whether the discount factor should be extrapolated for dates outside of the term structure’s range.

The implied forward rates can be obtained in a similar way.

InterestRate forwardRate(const Date& d1, const Date& d2, const DayCounter& resultDayCounter, Compounding comp, Frequency freq = Annual, bool extrapolate = false) const; InterestRate forwardRate(const Date& d, const Period& p, const DayCounter& resultDayCounter, Compounding comp, Frequency freq = Annual, bool extrapolate = false) const; InterestRate forwardRate(Time t1, Time t2, Compounding comp, Frequency freq = Annual, bool extrapolate = false) const;

For the forward rate there are three methods available, each taking the time interval in a different form. It can be given by two dates, by a date and a time period, or by two floating point numbers representing the time intervals. The behaviour is equivalent to the `zeroRate`

methods. For the first two methods, a separate day counter can be specified for the resulting `InterestRate`

. For the last version, the day counter of the term structure is used.

The forward rate is evaluated by calculating a compound rate using the ratio of the two discount rates, *a=d _{1}/d_{2}*. The interest rate is then generated by calling

`InterestRate::impliedRate()`

with the compound rate.As promised above, I have to now talk about the actual implementation of the discount factor calculation. This should be handled in the method `discountImpl`

which is declared abstract.

virtual DiscountFactor discountImpl(Time) const = 0;

Any concrete implementation of `YieldTermStructure`

must implement this method and return the discount factor for a specific time. When discountImpl() is called, the conversions from dates to units of time has already been done. The `enabled`

flag that was passed to the `discount`

method has been checked. This means that `discountImpl`

must assume that the discount factor has to be extrapolated if the time does not fall within the time for which the term structure can return values, i.e. if it is larger than the value returned by `maxTime`

.

Follow the author

]]>class LastFixingQuote : public Quote, public Observer

The class inherits from **Quote** but also from **Observer**. LastFixingQuote has to react to changes in the index value. It does this by observing the Index objects that it wraps.

LastFixingQuote(const boost::shared_ptr<Index>& index);

The constructor takes a shared pointer to an Index object. Once the index is set it cannot be changed again.

Real value() const; bool isValid() const;

The **value()** method and **isValid()** method are inherited from Quote and are implemented. **value()** returns the latest available fixing of the index. **isValid()** returns true if the index contains any fixings and false if the index does not contain fixings at all.

Date referenceDate() const;

The reference date is the last date in the time series contained in the index. Thus it is the date of the latest available fixing.

void update();

The **update()** method is implemented from the **Observer** class. It takes no action except calling **notifyObservers()** of the **Observable**. This makes sure that any observers of the quote are notified when the underlying index changes.

const boost::shared_ptr<Index>& index() const

Access to the index is granted through the index() method.

]]>virtual Real value() const = 0; virtual bool isValid() const = 0;

**value()** returns the value of a quote and **isValid()** returns true if the quote holds a valid value. Both these methods have to be implemented by a concrete sub-class of **Quote**.

**Quote** inherits from **Observable**. This makes it possible for observers to be notified when the value of a quote changes.

The simplest implementation of the abstract Quote class is SimpleQuote. SimpleQuote does not depend on an underlying quantity but stores the value of the quote directly. The value of the simple quote can be passed in the constructor.

SimpleQuote(Real value = Null<Real>());

This will initialize the quote with a given value. Note how the default value is **Null<Real>()**. This value indicates that the quote does not hold a valid value.

Real value() const; bool isValid() const;

These two methods have been implemented to return the stored value and to indicate if a valid value has been stored. The method **isValid()** will return false if the value which is stored in the SimpleQuote equals to **Null<Real>()**.

Real setValue(Real value = Null<Real>()); void reset();

The value of the simple quote can be modified by these two methods. setValue() will update the quote with a new value, or erase the value of the quote is no argument is given. The method will return the difference between the old value and the new value. If the value has changed any observers of the quote will be notified (see Observer and Observable). The convenience method **reset()** will erase any value stored in the quote. It is equivalent to calling **setValue()** without arguments.

The **Instrument** constructor, Instrument(), does not take any arguments. Derived classes will change this to include the data needed to evaluate the NPV of the instrument. Three methods give access to the evaluated data.

Real NPV() const; Real errorEstimate() const; const Date& valuationDate() const;

**NPV()** returns the net present value, **errorEstimate()** returns an estimate of the error contained in the NPV, and **valuationDate()** returns the date at which the instrument is valued. These methods are defined as const, since they apparently only access the underlying data. Because Instrument is a LazyObject these methods do, however, call the **performCalculations()** method which updates the internal data if necessary. Depending on the type of instrument, there may be additional data calculated along with the core values of NPV and the error estimate. These can be retrieved using the generic result() method.

template <typename T> T result(const std::string& tag) const;

This method takes the name of the result as argument and returns the value of that result. The caller of the **result()** method must know the type of the data that is stored under the given name. For example **Swaption**, which indirectly inherits from Instrument, also calculates the vega of the swaption. Imagine you have a **Swaption** object called **swaption**, then the vega can be accessed in the following way.

Real vega = swaption.result<Real>("vega");

Internally the additional results are stored in a **STL** map containing the names of the results as keys. The values are stored in **boost::any** objects. This allows any kind of data to be stored in the map but places the responsibility on the user to know exactly what type of data is stored. The complete map can be accessed through **additionalResults()**.

const std::map<std::string,boost::any>& additionalResults() const;

This method just returns a constant reference to the internally kept map of additional results. The results contained in the map depend on the implementation of the specific instrument class or, alternatively, the implementation of the pricing engine if this is used.

virtual bool isExpired() const = 0;

The NPV of an instrument is only calculated if the instrument has not expired. In the case that the instrument has expired the NPV and the error estimate are set to zero. All other additional results are not set to any default value but are deleted. **This means that any additional results of an instrument should not be accessed if it has expired!**

Because Instrument is derived from **LazyObject** it implements the method **performCalculations()**.

protected: virtual void performCalculations() const;

The default implementation of **performCalculations()** uses a pricing engine, derived from **PricingEngine**, to carry out the calculations. This default implementation makes sure that the right arguments are passed to the pricing engine and that they are validated for consistency before calculations are carried out. Finally it copies the results from the pricing engine to instrument’s internal fields. I will go into the details of pricing engines in a separate post. If the instrument should not use a pricing engine to calculate the results the **performCalculations()** method should be overridden. If a pricing engine is used it can be set with the **setPricingEngine()** method.

void setPricingEngine(const boost::shared_ptr<PricingEngine>&);

The following two methods are needed when a pricing engine is used. They are defined public but usually don’t need to be called because they are called from within **performCalculations()**.

virtual void setupArguments(PricingEngine::arguments*) const; virtual void fetchResults(const PricingEngine::results*) const;

The method setupArguments() sets up the arguments for the pricing engine to calculate the results. But note that the default implementation will throw an exception if called. **This means that setupArguments() must be implemented if a pricing engine is used!** If, on the other hand, the **performCalculations()** method is overridden and doesn’t use a pricing engine, then **setupArguments()** is not needed and doesn’t have to be implemented. The method **fetchResults()** fetches the NPV, the error estimate and any additional results from the pricing engine and stores them in the internal fields of the **Instrument** class.

protected: virtual void calculate() const;

Instrument overrides the **calculate()** method of the **LazyObject** class. As discussed in the post about lazy evaluation, this is not generally needed for lazy objects. **Instrument** modifes the existing behaviour to check for expiration of the instrument. If the instrument has expired then the NPV and error estimate is simply set to zero by calling and no evaluation takes place. If the instrument is not expired, the **calculate()** method of the parent class **LazyObject** is called to update the data.

protected: virtual void setupExpired() const;

This method is called when the instrument is expired to set the NPV and the error estimate to zero in case the instrument has expired. All the additional results are cleared from the internal storage and should not be accessed once this method has been called. If any other data should be cleared or any other action should be taken the setupExpired() method needs to be overridden. In that case it should be ensured, however, that the parent’s setupExpired() method is called.

]]>InterestRateIndex( const std::string& familyName, const Period& tenor, Natural fixingDays, const Currency& currency, const Calendar& fixingCalendar, const DayCounter& dayCounter);

The first parameter is the family name of the index. This could be something like **GBPLibor** or **Euribor**. The second parameter specifies the tenor of the particular index. The third parameter represents the fixing days, that is the number of business days between the fixing date and the value date. This is usually fixed for a particular family of interest rate indices, as is the currency of the index, the fixing calendar and the day counter. These last four parameters are usually specified in the constructor of a derived class and are not left for the user to specify.

std::string name() const;

This method overrides the **name()** method of the **Index** class. The name of the index is constructed from the information given in the constructor. It starts with the family name. For indices with a tenor of one day the ending depends on the settlement days and is either ON, TN or SN when **settlementDays** is 1, 2 or 3. For other tenors the tenor period is appended to the string in short format, e.g. 3M for a three month tenor.

Calendar fixingCalendar() const; bool isValidFixingDate(const Date& fixingDate) const;

The method **fixingCalendar()** will return the fixing calendar passed in the constructor. The method **isValidFixingDate()** is a shortcut and is equivalent to **fixingCalendar().isBusinessDay()**. Here are a few more inspectors.

std::string familyName() const; Period tenor() const; Natural fixingDays() const; const Currency& currency() const; const DayCounter& dayCounter() const;

All in all these return pretty much what you would expect. The values returned are just those that were passed in the constructor. The following two methods are used to convert between a value date and the corresponding fixing date.

Date fixingDate(const Date& valueDate) const; virtual Date valueDate(const Date& fixingDate) const;

The method **fixingDate()** calculates the fixing date from the value date. The value date lies a number of days after the fixing date, as specified by the **fixingDays** variable. To find the fixing date, the calendar is used to move backward by that number of days. Conversely, **valueDate()** will use the calendar to move forward by the number of business days as specified by **fixingDays** in order to find the value date from the fixing date. Note, how the **valueDate()** method is declared virtual whereas **fixingDate()** is not. This allows **valueDate()** to be re-defined in order to allow for specific conventions.

virtual Date maturityDate(const Date& valueDate) const = 0;

This method should calculate the maturity date from the value date of the index. Usually this will advance the value date by the tenor period. However, since there is a choice of business convention and other convention which influence the outcome, there is no default way of finding the maturity date. For this reason the method has been declared abstract and must be implemented in a concrete interest rate index class.

virtual Rate forecastFixing(const Date& fixingDate) const = 0; Rate pastFixing(const Date& fixingDate) const;

The method **forecastFixing()** is an abstract method that needs to be implemented by a concrete class in order to forecast fixings in the future. The method **pastFixing()**, on the other hand retrieves a previous fixing from the ones that have already been stored in the index.

Rate fixing(const Date& fixingDate, bool forecastTodaysFixing = false) const;

This method implements the method from Index and returns a fixing for a given date. If the date lies in the past then **pastFixing()** is called, if the date lies in the future then **forecastFixing()** is called. If fixing is called with today’s date (the evaluation date set in the global Settings) then the flag **forecastTodaysFixing** determines whether to forecast today’s fixing or not. If **forecastTodaysFixing** is true, then **fixing()** will always return a forecast value for today’s fixing. If, on the other hand, it is false then a past fixing is looked for first. If no past fixing for today is found the behavior depends on the result of

Settings::instance().enforcesTodaysHistoricFixings())

If this method indicates that today’s fixing should always be historic then an exception is thrown if it doesn’t exist. Otherwise **fixing()** will fall back to forecasting today’s fixing.