An Extensible C++ Framework for Stochastic Differential Equations

By Dave Burke

Engineers and scientists frequently use differential equations to model the effects of change, motion, and growth. For example, differential equations predict the motion of particles, changes in animal populations, and the movement of resources in a financial market. As models become more accurate and sophisticated, random or stochastic effects become more significant. In fact, the randomness in the system is sometimes the most interesting part of the model. To address such problems, a new area of mathematics called "stochastic calculus" has emerged. Stochastic calculus is actually a stochastic generalization of classical differential calculus. Just as you may use ordinary deterministic differential equations (ODEs) to model an ideal system with no noise, you use stochastic differential equations (SDEs) to model a more realistic system that includes the effects of noise. An SDE is a generalization of an ODE — when the noise is negligible, the equation reduces to an ODE. You can see this effect clearly in the conventional form of the SDE, which consists of two components:

• The drift term, the deterministic part.
• The diffusion term, the stochastic part.

Scientists who model real-world processes often begin with the assumption that noise is a passive, low-intensity background disturbance that may be averaged out to yield deterministic equations for the system. However, this practice of averaging noise effects often obscures distinctive and relevant behavior (for example, when the noise actually interacts with the system). Therefore, a stochastic model is the only realistic model for certain systems. A classic example of a case where noise interacts with the system is found in the problem of modeling the price of a stock option. Indeed, this application of SDEs initiated the entire options- and risk-management industry in finance, pioneered in the Nobel-prize-winning work of R. Merton and F. Black and M. Scholes [1][2]. By applying their work, you can calculate a fair price for a stock option at a given time — something that would not be possible using a set of deterministic equations that ignore stochastic effects. While it is likely that similar techniques will become pervasive in other commercial areas (insurance and marketing, for instance), there are other applications for SDEs in science and engineering, including radio astronomy, satellite-orbit stability, and genetics.

Numerical solvers are effective tools in an engineer's or scientist's armory for solving deterministic differential equations. The algorithms for implementing these solvers are well-known (Runge-Kutta, for instance) and widely available in software. Unfortunately, conventional solvers do not work for SDEs. In this article, I provide an introduction to stochastic equations, then present an extensible framework — written in C++ — for numerically solving SDEs. Given the accelerating interest in SDEs, the ability to extend this toolbox to accommodate new algorithms as they are developed is particularly attractive. The toolbox I describe here is also useful for existing applications that employ deterministic numerical methods and require the addition of the effects of noise.

The source code for the SDE toolbox (available at http://www.cuj.com/code/) includes three solvers of increasing accuracy and sophistication: Euler-Maruyama, Milstein, and Runge-Kutta. For brevity, I discuss only the Euler-Maruyama method in detail. To demonstrate the toolbox, I solve an SDE used in the Black-Scholes model for pricing stock options. The core classes for the toolbox are written in ANSI C++ for portability and have been tested with Microsoft Visual Studio C++ 6.0, Linux GCC 2.91, and 3.04. The toolbox is further packaged for Linux in "GNU style," facilitating a suite of command-line arguments and the ability to load dynamic systems at run time from shared libraries. This design reduces the need to recompile the main binary when changes are required.

The Basics

When I refer to stochastic differential equations, I mean equations of the form of Example 1(a), where a, which represents deterministic effects, is the drift component, and b, which represents stochastic effects, is the diffusion component. The term "noise" refers to Gaussian white noise and is replaced formally by the derivative of the so-called "Wiener" stochastic process, dW/dt, to give the usual SDE notation of Example 1(b).

The Wiener process W(t) is a random variable that changes over time and is also known as "Brownian Motion." For the one-dimensional case, the Wiener process describes the path taken by "a drunken ant with a clock and a coin": With each tick of the clock the ant tosses a coin and decides to take one step either left or right depending on the outcome. The result is a random walk, which on average goes nowhere (the expectation written as E(W) is 0) but has the potential to go further as time progresses (the variance written as E(W²) is equal to the elapsed time t, assuming the ant starts from 0). This concept of a path that "on average goes nowhere" but nevertheless has the potential to assume a direction due to the interaction of probabilistic effects is analogous to the "noise" within a stochastic system. Example 1(b) is is really just a symbolic way of representing the integral equation in Example 1(c).

The first integral in Example 1(c) is a normal (Riemann) integral that may be approximated by Example 1(d). In a similar fashion, the second integral may be interpreted as Example 1(e). This form is known as an "Ito" integral. What defines it as an Ito is the fact that we selected b at t_j (this is subtle because, for the normal integral, b can actually be selected anywhere between t_j+1 and t_j). The Ito integral is part of what is known as the "Ito calculus" and includes an alternative to the usual chain rule known as "Ito's formula." The trick to numerically solving Ito SDEs is representing the Wiener increment W(t_j+1)-W(t_j) algorithmically. As with deterministic differential equations, it is typically straightforward to find explicit solutions to linear equations, where (in the stochastic case) terms in W(t) are present. However, and again as in the deterministic case, many equations do not possess explicit solutions (for example, most nonlinear equations); hence, you often must resort to numerical techniques to seek solutions. For the purposes of the toolbox, vector SDEs are employed so that X, W, and a are d-dimensional vectors, and b is a d×d dimensional matrix. This approach lets you solve multidimensional systems.

Architecture

Figure 1 illustrates the object-oriented architecture of the toolbox. Solver implementations are derived from the base class CSolver and instantiated via the Parameterized Factory CSolverFactory [4]. The CSolverSettings class aggregates solver-independent settings. To allow solvers to be developed independently of the dynamic system, you are solving (and vice versa), a generic interface called CDynSystem represents dynamic systems. This class has the method GetDriftAndDiffusion(), which returns the a and b components. This design is noteworthy as it differs from traditional (deterministic) solvers, which often employ a function to return the entire right hand side of the equation you're solving. Other algorithms for solving SDEs manipulate the drift and diffusion components differently; hence, they need to be obtained from the dynamic system separately.

To ease the development of new solvers, the Template Method pattern [4] is employed so that CSolver defines a skeleton of the solution procedure, deferring the specific details of advancing the solution to subclasses. Figure 2 is a sequence diagram of the solution procedure.

A solver is invoked by calling the Run() method. This calls Initialise(), which seeds the random-number generator (either by user-supplied values or values derived from the system time), creates the Wiener process realization, and initializes memory for the solution. If specified, the Initialise() method writes the Wiener process realization to disk as well. For each time step, the drift and diffusion components of the system to be solved are retrieved via the GetDriftAndDiffusion() and the Advance() methods of the particular solver. The solution process terminates with a call to the Finish() method, which writes the solution to disk; see Listing 1.

Euler-Maruyama Solver

The Euler-Maruyama method is probably the simplest solver for SDEs and is an extension of the Euler method used for deterministic differential equations. The scheme follows from Example 1(c) if you consider a uniform discretization of N equal time steps from time 0 to time T. This yields a time step of =_n+1-_n=T/N. The scheme may be written as Example 2(a), where W_n=W(_n+1)-W(_n). For b=0, you obtain the usual deterministic Euler scheme. The random variables in W_n=W(_n+1)-W(_n) may be obtained from a Gaussian distribution of mean 0 and variance t. The CSolver base class creates these random variables during the Initialise() method by invoking an efficient Box-Muller transformation of a uniformly generated sequence of random numbers [6]. The latter are generated via the "Mersenne Twister" [5], a high-quality random-number generator. The random-number generator in the Standard C Library generally performs poorly when employed in stochastic modeling settings. To simplify the implementation and prevent the introduction of errors, the Wiener process is calculated at m_nStepSize increments, and the solver computes the solution at increments of m_nStepSize*m_nMultiplier, thus avoiding the need for interpolation.

To compare the quality of different solvers, it is instructive to introduce a measure of accuracy. The strong order of convergence is calculated from Example 2(b) for all step sizes (0,1) for some constant K, and is a useful consideration in applications that require sample path information (such as for the generation of a stock market price scenario). In some applications, only statistical information about the solution, such as the mean and standard deviation, may be required, and for those cases, a weaker form of convergence measure is more relevant [3]. The Euler-Maruyama method has a strong order of convergence of 0.5. The method is obtained from a truncation of a stochastic version of the Taylor series, and by including higher order terms you obtain higher orders of convergence — the Milstein solver has strong order 1.0 and the Runge-Kutta solver has strong order 1.5.

An Application

To demonstrate an application of the toolbox, I numerically solve an SDE from the Black-Scholes model for pricing derivatives (see the accompanying text box entitled "The Black-Scholes Model"). A core assumption of the model is that a stock price may be described by the SDE in Example 3(a), where m represents the average rate of growth of the stock, and s is the uncertainty or volatility of the stock. This equation is also known as "Geometric Brownian Motion." Example 3(a) is a straightforward SDE that possesses the explicit solution given by Example 3(b). This explicit solution lets you compare the actual solution of Example 3(a) with a numerical solution using different solvers from the toolbox. The toolbox has a setting that enables the W_t series to be written to disk; hence you can easily plot Example 3(b). Figure 3 illustrates Example 3(b) with the corresponding Euler-Maruyama solution of Example 3(a). Figure 4 is the same as Figure 3, except I used the Runge-Kutta solver, which clearly produces a more accurate solution for the same step size (in fact, the solutions coincide perfectly). Finally, Figure 5 illustrates three simulations of Example 3(a) using a different realization of the Wiener process each time and an average of 100 simulations exposing the underlying growth.

Conclusions

The toolbox presented here provides implementations for the Euler-Maruyama (strong order 0.5), Milstein (strong order 1.0), and Runge-Kutta (strong order 1.5) solvers, along with the ability to quickly extend the framework for new solvers (subject to the constraint that they operate explicitly and do not require derivatives). I used the example of a simple (linear) SDE, but the real power of these numerical techniques comes with their application to complex nonlinear equations. I've used this toolbox to study prescribed stochastic limit cycle oscillators (see CEllipseDynSystem in the included source for example). Hopefully, this article has given you a small taste of stochastic differential equations.

References

[1] Black, F., and M. Scholes, "The Pricing of Options and Corporate Liabilities," J. Political Economy, 1973, p. 637-659.

[2] Merton, R.,"The Theory of Rational Option Pricing," Bell Journal of Economics and Management Science, 1973, p. 141-183.

[3] Kloeden, P.E., and E. Platen. Numerical Solution of Stochastic Differential Equations, Third Edition, Springer-Verlag, 1999.

[4] Gamma, E., R. Helm, R. Johnson, and J. Vlissides. Design Patterns: Elements of Reusable Object-Oriented Software, Addison-Wesley, 1995.

[5] Matsumoto, M., and T. Nishimura, "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator," ACM Transactions on Modeling and Computer Simulation, 1998, p. 3.

[6] Box, G.E.P, and M.E. Muller, "A Note on the Generation of Random Normal Deviates," Annals Math. Stat, 1958, p. 610-611.