# Chapter 6, part 2 of n: stochastic processes

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 *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 buttons for that are in the footer. 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.

Listing 6.4: Partial implementation of the `StochasticProcess`

class.

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 `Array`

s 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 `Real`

s into and from `Array`

s, 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.

#### Bibliography

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