Chapter 6, part 2 of n: stochastic processes
Share:
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.
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.
Share: