Solution of Stochastic Differential Equations in Finance

Table 1. Average error at T = 8 of approximate solutions of (4). The error scales as ∆t^1/2 for Euler-Maruyama and ∆t for Milstein.

∆t	Euler-Maruyama	Milstein
2−1	0.169369	0.063864
2−2	0.136665	0.035890
2−3	0.086185	0.017960
2−4	0.060615	0.008360
2−5	0.048823	0.004158
2−6	0.035690	0.002058
2−7	0.024277	0.000981
2−8	0.016399	0.000471
2−9	0.011897	0.000242
2−10	0.007913	0.000122

10⁰
₁₀−1

mean error
₁₀−2
₁₀−3
₁₀−4
10⁻⁴ 10⁻² 10⁰
stepsize t
Fig. 3. Error in the Euler-Maruyama and Milstein methods. Solution paths are computed for the geometric Brownian motion equation (15) and are compared to the correct X(T ) given by (5). The absolute difference is plotted versus stepsize h for the two different methods. The Euler-Maruyama errors are plotted as circles and the Milstein error as squares. Note the slopes 1/2 and 1, respectively, on the log-log plot.
In the stochastic differential equation context, the same trade can be made with the Milstein method, resulting in a strong order 1 method that requires evaluation of b(X) at two places on each step. A heuristic derivation can be carried out by making the replacement
b(w_i + b(w_i)^√∆t_i ) − b(w_i)

b_x(w_i) ≈
b(w )^√∆t

i i
in the Milstein formula (14), which leads to the following method (Rumelin, 1982):

Strong Order 1.0 Runge-Kutta Method

w₀ = X₀
w_i+1 = w_i + a(w_i)∆t_i + b(w_i)∆W_i
2 [b(w_i + b(w_i) ∆t_i) − b(w_i)](∆W_i − ∆t_i)/ ∆t_i
+ 1 2

The orders of the methods introduced here for SDEs, 1/2 for Euler- Maruyama and 1 for Milstein and the Runge-Kutta counterpart, would be considered low by ODE standards. Higher-order methods can be developed for SDEs, but become much more complicated as the order grows. As an ex- ample, consider the strong order 1.5 scheme for the SDE (13) proposed in Platen and Wagner (1982):

Strong Order 1.5 Taylor Method

w₀ = X₀
w_i+1 = w_i + a∆t_i + b∆W_i + ¹ bb_x(∆W ² − ∆t_i)
2 ⁱ
+ a_y σ∆Z_i + ¹ (aa_x + ¹ b²a_xx)∆t²
2 2 ⁱ

2
+ (ab_x + ¹ b²b_xx)(∆W_i∆t_i − ∆Z_i)
₊ 1 2 1 2
₂ b(bb_xx + b_x)( ₃ ∆W_i − ∆t_i)∆W_i (17)

i
where partial derivatives are denoted by subscripts, and where the additional random variable ∆Z_i is normally distributed with mean 0, variance E(∆Z²) =
¹ 3 ¹ 2
₃ ∆t_i and correlated with ∆W_i with covariance E(∆Z_i∆W_i) = ₂ ∆t_i . Note that ∆Z_i can be generated as

2
∆Z_i = ¹ ∆t_i(∆W_i + ∆V_i/^√3) where ∆V_i is chosen independently from ^√∆t_iN (0, 1).
Whether higher-order methods are needed in a given application depends
on how the resulting approximate solutions are to be used. In the ordinary differential equation case, the usual assumption is that the initial condition and the equation are known with accuracy. Then it makes sense to calculate the solution as closely as possible to the same accuracy, and higher-order methods are called for. In the context of stochastic differential equations, in particular if the initial conditions are chosen from a probability distribution, the advantages of higher-order solvers are often less compelling, and if they come with added computational expense, may not be warranted.

4. Weak convergence of SDE solvers

Strong convergence allows accurate approximations to be computed on an individual realization basis. For some applications, such detailed pathwise

information is required. In other cases, the goal is to ascertain the probability distribution of the solution X(T ), and single realizations are not of primary interest.

Weak solvers seek to fill this need. They can be simpler than corresponding strong methods, since their goal is to replicate the probability distribution only. The following additional definition is useful.
A discrete-time approximation w_∆t with step-size ∆t is said to converge weakly to the solution X(T ) if

lim
∆t→0
E{f (w_∆t(T ))} = E{f (X(T ))}

for all polynomials f (x). According to this definition, all moments converge as ∆t → 0. If the stochastic part of the equation is zero and the initial value is deterministic, the definition agrees with the strong convergence definition, and the usual ordinary differential equation definition.
Weakly convergent methods can also be assigned an order of convergence. We say that a the solver converges weakly with order m if the error in the moments is of mth order in the stepsize, or
|E{f (X(T ))} − E{f (w_∆t(T ))}| = O((∆t)^m)
for sufficiently small stepsize ∆t.
In general, the rates of weak and strong convergence do not agree. Unlike the case of ordinary differential equations, where the Euler method has order 1, the Euler-Maruyama method for SDEs has strong order m = 1/2. However, Euler-Maruyama is guaranteed to converge weakly with order 1.
Higher order weak methods can be much simpler than corresponding strong methods, and are available in several different forms. The most direct approach is to exploit the Ito-Taylor expansion (Kloeden and Platen, 1992), the Ito calculus analogue of the Taylor expansion of deterministic functions. An example SDE solver that converges weakly with order 2 is the following: