Hello everybody.

This post is the second in a series that covers chapter 6 of my book, that is, the Monte Carlo framework in QuantLib. The first post is here. In this post: stochastic processes.

Also, as I mentioned last week, you can now register for my next

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.

This post is the second in a series that covers chapter 6 of my book, that is, the Monte Carlo framework in QuantLib. The first post is here. In this post: stochastic processes.

Also, as I mentioned last week, you can now register for my next

*Introduction to QuantLib Development*course (which is going to be in London, September 22nd to 24th). You can get an early-bird discount until the end of May.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.

### Stochastic processes

The whole point of using RNGs is, of course, to use the generated random numbers to drive a stochastic process. The interface of the StochasticProcess class (displayed in listing 6.4) was designed with this application in mind—and it shows.

$$

\mathrm{d}\mathbf{x} = \mathbf{\mu}(t, \mathbf{x}) \mathrm{d}t

+ \mathbf{\sigma}(t, \mathbf{x}) \cdot \mathrm{d}\mathbf{w}

$$

However, the interface of the class deals more with the discretization of the process over the time steps of a sampled path.

Straight away, the class defines a polymorphic inner class discretization. It is an instance of the Strategy pattern, whose purpose is to make it possible to change the way a process is discretely sampled. Its methods take a StochasticProcess instance, a starting point \( (t,\mathbf{x}) \), and a time step \( \Delta t \) and returns a discretization of the process drift and diffusion. (By this time, you'll have surely noticed that there's no provision for jumps. Yes, I know. Hold that thought; I'll have more to say on this later.)

Back to the StochasticProcess interface. The class defines a few inspectors. The size method returns the number of variables modeled by the process. For instance, it would return 1 for a regular Black-Scholes process with a single underlying, or 2 for a stochastic-volatility model, or N for the joint process of N correlating underlying assets each following a Black-Scholes process. The factors method returns the number of random variates used to drive the variables; it has a default implementation returning the same number as size, although on afterthought I'm not sure that it was a great idea; we might have required to specify that explicitly. To round up the inspectors, the initialValues method returns the values of the underlying variables at t=0; in other words, the present values, in which no randomness is involved. As in most StochasticProcess methods, the result is wrapped in a Disposable class template. This is a way to avoid a few copies when returning arrays or matrices; if you're interested in the details, you can have a look at appendix A. When we can rely on compilers supporting the C++11 standard, we'll be able to replace those with r-value references [1]; but as usual, it will be a few years before we can do this (where "this" also includes being able to drop support for older compilers with a clear conscience, and without leaving users outside in the rain).

The next pair of methods (also inspectors, in a way) return more detailed information on the process. The drift method return the $\mu(t,\mathbf{x})$ term in the previous formula, while the diffusion method returns the $\sigma(t,\mathbf{x})$ term. Obviously, for dimensions to match, the drift term must be a N-dimensional array and the diffusion term a M × N matrix, with M and N being the dimensions returned by the factors and size methods, respectively. Also, the diffusion matrix must include correlation information.

The last set of methods deals with the discretization of the process over a path. The expectation method returns the expectation values of the process variables at time \( t+\Delta t \), assuming they had values \( \mathbf{x} \) at time \( t \). Its default implementation asks the stored discretization instance for the integrated drift of the variables over \( \Delta t \) and applies it to the starting values (apply is another virtual method that takes a value \( \mathbf{x} \) and a variation \( \Delta\mathbf{x} \) and unsurprisingly returns \( \mathbf{x} + \Delta\mathbf{x} \). Why we needed a polymorphic method to do this, you ask? It's a sad story that will be told in a short while). The stdDeviation method returns the diffusion term integrated over \( t+\Delta t \) and starting at \( (t, \mathbf{x}) \), while the covariance method returns its square; in their default implementation, they both delegate the calculation to the discretization instance. Finally, the evolve method provides some kind of higher-level interface; it takes a starting point \( (t, \mathbf{x}) \), a time step \( \Delta t \), and an array \( \mathbf{w} \) of random Gaussian variates, and returns the end point \( \mathbf{x'} \). Its default implementation asks for the standard-deviation matrix, multiplies it by the random variates to yield the random part of the simulated path, and applies the resulting variations to the expectation values of the variables at \( t+\Delta t \) (which includes the deterministic part of the step).

As you've seen, the methods in this last set are all given a default implementation, sometimes calling methods like drift or diffusion by way of the discretization strategy and sometimes calling directly other methods in the same set; it's a kind of multiple-level Template Method pattern, as it were. This makes it possible for classes inheriting from StochasticProcess to specify their behavior at different levels.

At the very least, such a class must implement the drift and diffusion methods, which are purely virtual in the base class. This specifies the process at the innermost level, and is enough to have a working stochastic process; the other methods would use their default implementation and delegate the needed calculations to the contained discretization instance. The derived class can prescribe a particular discretization, suggest one by setting it as a default, or leave the choice entirely to the user.

At a somewhat outer level, the process class can also override the expectation and stdDeviation methods; this might be done, for instance, when it's possible to write an closed formula for their results. In this case, the developer of the class can decide either to prevent the constructor from taking a discretization instance, thus forcing the use of the formula; or to allow the constructor to take one, in which case the overriding method should check for its presence and possibly forward to the base-class implementation instead of using its own. The tasks of factoring in the random variates and of composing the integrated drift and diffusion terms are still left to the evolve method.

At the outermost level, the derived class can override the evolve method and do whatever it needs—even disregard the other methods in the StochasticProcess interface and rely instead on ones declared in the derived class itself. This can serve as a Hail Mary pass for those process which don't fit the interface as currently defined. For instance, a developer wanting to add stochastic jumps to a process might add them into evolve, provided that a random-sequence generator can generate the right mixture of Gaussian and Poisson variates to pass to the method.

Wouldn't it be better to support jumps in the interface? Well, yes, in the sense that it would be better if a developer could work

On the one hand, I'm not sure what the right interface should be. For instance, can we settle on a single generalized formula including jumps? I wouldn't want to release an interface just to find out later that it's not the correct one; which is likely, given that at this time we don't have enough code to generalize from.

On the other hand, what happens when we try to integrate the next feature? Do we add it to the interface, too? Or should we try instead to generalize by not specializing; that is, by declaring fewer methods in the base-class interface, instead of more? For instance, we could keep evolve or something like it, while moving drift, diffusion and their likes in subclasses instead (by the way, this might be a way to work on a jump interface without committing too early to it: writing some experimental classes that define the new methods and override evolve to call them. Trying to make them work should provide some useful feedback).

But I went on enough already. In short: yes, Virginia, the lack of jumps in the interface is a bad thing. However, I think we'd be likely to make a poor job if we were to add them now. I, for one, am going to stay at the window and see if anybody comes along with some working code. If you wanted to contribute, those experimental classes would be most appreciated.

One last detour before we tackle a few examples. When implementing a one-dimensional process, the Arrays and their likes tend to get in the way; they make the code more tiresome to write and less clear to read. Ideally, we'd want an interface like the one declared by StochasticProcess but with methods that take and return simple real numbers instead of arrays and matrices. However, we can't declare such methods as overloads in StochasticProcess, since they don't apply to a generic process; and neither we wanted to have a separate hierarchy for 1-D processes, since they belong with all the others.

The library solves his problem by providing the StochasticProcess1D class, shown in listing 6.5. On the one hand, this class declares an interface that mirrors exactly the StochasticProcess interface, even down to the default implementations, but uses the desired Real type. On the other hand, it inherits from StochasticProcess (establishing that a 1-D process "is-a" stochastic process, as we wanted) and implements its interface in terms of the new one by boxing and unboxing Reals into and from Arrays, so that the developer of a 1-D process needs not care about it: in short, a straightforward application of the Adapter pattern. Because of the 1:1 correspondence between scalar and vectorial methods (as can be seen in the listing) the 1-dimensional process will work the same way when accessed from either interface; also, it will have the same possibilities of customization provided by the StochasticProcess class and described earlier in this section.

Listing 6.4: Partial implementation of the StochasticProcess class.

```
class StochasticProcess : public Observer,
public Observable {
public:
class discretization {
public:
virtual ~discretization() {}
virtual Disposable<Array> drift(const StochasticProcess&,
Time t0, const Array& x0,
Time dt) const = 0;
// same for diffusion etc.
};
virtual ~StochasticProcess() {}
virtual Size size() const = 0;
virtual Size factors() const {
return size();
}
virtual Disposable<Array> initialValues() const = 0;
virtual Disposable<Array> drift(Time t,
const Array& x) const = 0;
virtual Disposable<Matrix> diffusion(Time t,
const Array& x)
const = 0;
virtual Disposable<Array> expectation(Time t0,
const Array& x0,
Time dt) const {
return apply(x0, discretization_->drift(*this,t0,x0,dt));
}
virtual Disposable<Matrix> stdDeviation(Time t0,
const Array& x0,
Time dt) const {
return discretization_->diffusion(*this,t0,x0,dt);
}
virtual Disposable<Matrix> covariance(Time t0,
const Array& x0,
Time dt) const {
return discretization_->covariance(*this,t0,x0,dt);
}
virtual Disposable<Array> evolve(Time t0,
const Array& x0,
Time dt,
const Array& dw) const {
return apply(expectation(t0,x0,dt),
stdDeviation(t0,x0,dt)*dw);
}
virtual Disposable<Array> apply(const Array& x0,
const Array& dx) const {
return x0 + dx;
}
};
```

The concept modeled by the class is a somewhat generic n-dimensional stochastic process described by the equation$$

\mathrm{d}\mathbf{x} = \mathbf{\mu}(t, \mathbf{x}) \mathrm{d}t

+ \mathbf{\sigma}(t, \mathbf{x}) \cdot \mathrm{d}\mathbf{w}

$$

However, the interface of the class deals more with the discretization of the process over the time steps of a sampled path.

Straight away, the class defines a polymorphic inner class discretization. It is an instance of the Strategy pattern, whose purpose is to make it possible to change the way a process is discretely sampled. Its methods take a StochasticProcess instance, a starting point \( (t,\mathbf{x}) \), and a time step \( \Delta t \) and returns a discretization of the process drift and diffusion. (By this time, you'll have surely noticed that there's no provision for jumps. Yes, I know. Hold that thought; I'll have more to say on this later.)

Back to the StochasticProcess interface. The class defines a few inspectors. The size method returns the number of variables modeled by the process. For instance, it would return 1 for a regular Black-Scholes process with a single underlying, or 2 for a stochastic-volatility model, or N for the joint process of N correlating underlying assets each following a Black-Scholes process. The factors method returns the number of random variates used to drive the variables; it has a default implementation returning the same number as size, although on afterthought I'm not sure that it was a great idea; we might have required to specify that explicitly. To round up the inspectors, the initialValues method returns the values of the underlying variables at t=0; in other words, the present values, in which no randomness is involved. As in most StochasticProcess methods, the result is wrapped in a Disposable class template. This is a way to avoid a few copies when returning arrays or matrices; if you're interested in the details, you can have a look at appendix A. When we can rely on compilers supporting the C++11 standard, we'll be able to replace those with r-value references [1]; but as usual, it will be a few years before we can do this (where "this" also includes being able to drop support for older compilers with a clear conscience, and without leaving users outside in the rain).

The next pair of methods (also inspectors, in a way) return more detailed information on the process. The drift method return the $\mu(t,\mathbf{x})$ term in the previous formula, while the diffusion method returns the $\sigma(t,\mathbf{x})$ term. Obviously, for dimensions to match, the drift term must be a N-dimensional array and the diffusion term a M × N matrix, with M and N being the dimensions returned by the factors and size methods, respectively. Also, the diffusion matrix must include correlation information.

The last set of methods deals with the discretization of the process over a path. The expectation method returns the expectation values of the process variables at time \( t+\Delta t \), assuming they had values \( \mathbf{x} \) at time \( t \). Its default implementation asks the stored discretization instance for the integrated drift of the variables over \( \Delta t \) and applies it to the starting values (apply is another virtual method that takes a value \( \mathbf{x} \) and a variation \( \Delta\mathbf{x} \) and unsurprisingly returns \( \mathbf{x} + \Delta\mathbf{x} \). Why we needed a polymorphic method to do this, you ask? It's a sad story that will be told in a short while). The stdDeviation method returns the diffusion term integrated over \( t+\Delta t \) and starting at \( (t, \mathbf{x}) \), while the covariance method returns its square; in their default implementation, they both delegate the calculation to the discretization instance. Finally, the evolve method provides some kind of higher-level interface; it takes a starting point \( (t, \mathbf{x}) \), a time step \( \Delta t \), and an array \( \mathbf{w} \) of random Gaussian variates, and returns the end point \( \mathbf{x'} \). Its default implementation asks for the standard-deviation matrix, multiplies it by the random variates to yield the random part of the simulated path, and applies the resulting variations to the expectation values of the variables at \( t+\Delta t \) (which includes the deterministic part of the step).

As you've seen, the methods in this last set are all given a default implementation, sometimes calling methods like drift or diffusion by way of the discretization strategy and sometimes calling directly other methods in the same set; it's a kind of multiple-level Template Method pattern, as it were. This makes it possible for classes inheriting from StochasticProcess to specify their behavior at different levels.

At the very least, such a class must implement the drift and diffusion methods, which are purely virtual in the base class. This specifies the process at the innermost level, and is enough to have a working stochastic process; the other methods would use their default implementation and delegate the needed calculations to the contained discretization instance. The derived class can prescribe a particular discretization, suggest one by setting it as a default, or leave the choice entirely to the user.

At a somewhat outer level, the process class can also override the expectation and stdDeviation methods; this might be done, for instance, when it's possible to write an closed formula for their results. In this case, the developer of the class can decide either to prevent the constructor from taking a discretization instance, thus forcing the use of the formula; or to allow the constructor to take one, in which case the overriding method should check for its presence and possibly forward to the base-class implementation instead of using its own. The tasks of factoring in the random variates and of composing the integrated drift and diffusion terms are still left to the evolve method.

At the outermost level, the derived class can override the evolve method and do whatever it needs—even disregard the other methods in the StochasticProcess interface and rely instead on ones declared in the derived class itself. This can serve as a Hail Mary pass for those process which don't fit the interface as currently defined. For instance, a developer wanting to add stochastic jumps to a process might add them into evolve, provided that a random-sequence generator can generate the right mixture of Gaussian and Poisson variates to pass to the method.

Wouldn't it be better to support jumps in the interface? Well, yes, in the sense that it would be better if a developer could work

*with*the interface rather than*around*it. However, there are two kinds of problems that I see about adding jumps to the StochasticProcess class at this time.On the one hand, I'm not sure what the right interface should be. For instance, can we settle on a single generalized formula including jumps? I wouldn't want to release an interface just to find out later that it's not the correct one; which is likely, given that at this time we don't have enough code to generalize from.

On the other hand, what happens when we try to integrate the next feature? Do we add it to the interface, too? Or should we try instead to generalize by not specializing; that is, by declaring fewer methods in the base-class interface, instead of more? For instance, we could keep evolve or something like it, while moving drift, diffusion and their likes in subclasses instead (by the way, this might be a way to work on a jump interface without committing too early to it: writing some experimental classes that define the new methods and override evolve to call them. Trying to make them work should provide some useful feedback).

But I went on enough already. In short: yes, Virginia, the lack of jumps in the interface is a bad thing. However, I think we'd be likely to make a poor job if we were to add them now. I, for one, am going to stay at the window and see if anybody comes along with some working code. If you wanted to contribute, those experimental classes would be most appreciated.

One last detour before we tackle a few examples. When implementing a one-dimensional process, the Arrays and their likes tend to get in the way; they make the code more tiresome to write and less clear to read. Ideally, we'd want an interface like the one declared by StochasticProcess but with methods that take and return simple real numbers instead of arrays and matrices. However, we can't declare such methods as overloads in StochasticProcess, since they don't apply to a generic process; and neither we wanted to have a separate hierarchy for 1-D processes, since they belong with all the others.

The library solves his problem by providing the StochasticProcess1D class, shown in listing 6.5. On the one hand, this class declares an interface that mirrors exactly the StochasticProcess interface, even down to the default implementations, but uses the desired Real type. On the other hand, it inherits from StochasticProcess (establishing that a 1-D process "is-a" stochastic process, as we wanted) and implements its interface in terms of the new one by boxing and unboxing Reals into and from Arrays, so that the developer of a 1-D process needs not care about it: in short, a straightforward application of the Adapter pattern. Because of the 1:1 correspondence between scalar and vectorial methods (as can be seen in the listing) the 1-dimensional process will work the same way when accessed from either interface; also, it will have the same possibilities of customization provided by the StochasticProcess class and described earlier in this section.

Listing 6.5: Partial implementation of the StochasticProcess1D class.

```
class StochasticProcess1D : public StochasticProcess {
public:
class discretization {
public:
virtual ~discretization() {}
virtual Real drift(const StochasticProcess1D&,
Time t0, Real x0, Time dt) const = 0;
// same for diffusion etc.
};
virtual Real x0() const = 0;
virtual Real drift(Time t, Real x) const = 0;
virtual Real diffusion(Time t, Real x) const = 0;
virtual Real expectation(Time t0, Real x0, Time dt) const {
return apply(x0, discretization_->drift(*this, t0,
x0, dt));
}
virtual Real stdDeviation(Time t0, Real x0, Time dt) const {
return discretization_->diffusion(*this, t0, x0, dt);
}
virtual Real variance(Time t0, Real x0, Time dt) const;
virtual Real evolve(Time t0, Real x0,
Time dt, Real dw) const {
return apply(expectation(t0,x0,dt),
stdDeviation(t0,x0,dt)*dw);
}
private:
Size size() const { return 1; }
Disposable<Array> initialValues() const {
return Array(1, x0());
}
Disposable<Array> drift(Time t, const Array& x) const {
return Array(1, drift(t, x[0]));
}
// ...the rest of the \texttt{StochasticProcess} interface...
Disposable<Array> evolve(Time t0, const Array& x0,
Time dt, const Array& dw) const {
return Array(1, evolve(t0,x0[0],dt,dw[0]));
}
};
```

#### Bibliography

[1] H.E. Hinnant, B. Stroustrup and B. Kozicki,*A Brief Introduction to Rvalue References*, C++ Standards Committee Paper N2027, 2006.