National Center for Theoretical Sciences (NCTS) at the National Tsing Hua University

STOCHASTIC DIFFERENTIAL EQUATIONS ARE BECOMING INCREASINGLY MORE POPULAR

(1) Random effects in physical, biological, and financial problems can be modeled using stochastic differential equations.

(2) Stochastic models are considered to be more realistic for many prob-lems.

(3) However, the development of accurate stochastic differential equation models is not well understood.

(4) In this series of lectures, stochastic differential equations are introduced, numerical methods are explained, and a procedure is described for deriving stochastic differential equations.

THESE LECTURES ARE DIVIDED INTO FOUR PARTS

Part I: Random variables and stochastic processes are reviewed and sto-chastic integrals are introduced.

Part II: The theory and approximation of Itˆo stochastic differential equa-tions (SDEs) are studied.

Part III: A procedure is explained for deriving accurate SDE models. SDEs are derived for several problems in biology, physics, and finance.

Part IV: The derivation procedure is extended to stochastic partial dif-ferential equations (SPDEs). SPDEs are derived for several problems in physics and biology.

STOCHASTIC DIFFERENTIAL EQUATIONS

It is proved that a unique solution to an SDE exists in Hilbert space H_SP. Several properties of stochastic differential equations are derived.

Exact solutions and moments of the solution are found.

Numerical methods for approximating solutions of SDEs are described.

The forward Kolmogorov partial differential equation is derived.

A procedure is described for estimating parameters in an SDE.

IT ˆO STOCHASTIC DIFFERENTIAL EQUATIONS We study

X(t, ω) = X(0, ω) + Z t

0

f (s, X(s, ω)) ds + Z t

0

g(s, X(s, ω)) dW (s, ω) for 0 ≤ t ≤ T where X(0, ·) ∈ H^RV or in differential form

dX(t, ω) = f (t, X(t, ω)) dt + g(t, X(t, ω)) dW (t, ω).

f is called the drift coefficient and g is called the diffusion coefficient.

It is assumed f and g satisfy:

Condition (c6): |f(t, x)−f(s, y)|² ≤ k(|t−s|+|x−y|²) for 0 ≤ s, t ≤ T and x, y ∈ R.

Condition (c7): |f(t, x)|² ≤ k(1 + |x|²) for 0 ≤ t ≤ T and x ∈ R.

PROPERTIES OF f AND g THAT SATISFY (c6) AND (c7)

Let u(t) = f (t, X(t)). Then, (c7) implies that u ∈ H^SP when X ∈ H^SP and, of course, the same holds for function g. Indeed, for X ∈ H^SP,

kuk²SP = Z T

0 E|f(t, X(t)|²dt ≤ Z T

0 k(1 + E|X(t)|²) dt ≤ kT + kkXk²SP.

Condition (c6) implies that kf(t¹, X) −f(t¹, Y )k^SP < ǫ when kX −Y k^SP < ǫ/k^1/2 and kf(t¹, X) − f(t², X)k^SP < ǫ when |t² − t¹| < ǫ²/kT .

To see the first inequality, consider kf(t¹, X) − f(t¹, Y )k²SP =

Z T

0 E|f(t¹, X(t)) − f(t¹, Y (t))|²dt

≤

Z T

0 kE|X(t) − Y (t)|²dt = kkX − Y k²SP.

Thus, if kX − Y k^SP < ǫ/k^1/2, then kf(t¹, X) − f(t¹, Y )k^SP < ǫ. To see the second inequality, consider

kf(t¹, X) − f(t², X)k²SP =

Z T

0 E|f(t¹, X(t)) − f(t², X(t))|²dt

≤ Z T

0 kE|t² − t¹| dt = kT |t² − t¹|.

Thus, if |t² − t¹| < ǫ²/kT then kf(t¹, X) − f(t², X)k^SP < ǫ.

EXISTENCE OF A UNIQUE SOLUTION

We now prove existence and uniqueness for X ∈ H^SP that solves the SDE when f and g satisfy conditions (c6) and (c7). To show existence and uniqueness, a Cauchy sequence of functions in HSP is constructed. The limit of this sequence is the solution.

Let X₀(t) = X(0), where X(0) ∈ H^RV is the given initial condition. Define This inequality follows from the following argument.

kX¹ − X⁰k²SP =

Continuing this procedure, the sequence {Xⁿ}^∞n=1 ⊂ H^SP is defined as

To see that this sequence is Cauchy in HSP, X_n+1(t) − Xⁿ(t) =

Let an(t) = E|Xⁿ⁺¹(t) − Xⁿ(t)|², and then, by the previous inequality, This latter inequality implies that

an(t) ≤ Lⁿtⁿ⁻¹

(n − 1)!kX¹ − X⁰k²SP

and it follows that

kXⁿ⁺¹ − Xⁿk²SP ≤ LⁿTⁿ

n! kX¹ − X⁰k²SP.

Consider m > n. By repeated application of the triangle inequality, kX^m − Xⁿk^SP ≤ kX^m − Xm−1k^SP + kXm−1 − Xm−2k^SP +· · · + kXⁿ⁺¹ − Xⁿk^SP

From the previous inequality, as kX¹ − X⁰k^SP is bounded, given ǫ > 0 there exists an N > 0 such that kXⁿ − X^mk^SP < ǫ when n, m > N . Hence, {Xⁿ}^∞n=1 is Cauchy in H_SP and converges to a unique X ∈ H^SP.

Now let Y ∈ H^SP where Y satisfies the relation Y (t) = −X(t) + X(0) +

Z t 0

f (s, X(s)) ds + Z t

0

g(s, X(s)) dW (s).

Because

0 = −Xⁿ(t) + X(0) + Z t

0

f (s, X_n−1(s)) ds + Z t

0

g(s, X_n−1(s)) dW (s) and kX − Xⁿk^SP → 0 as n → ∞, it is clear that kY k^SP = 0.

So X(t) is the unique solution in HSP.

PROPERTIES OF SOLUTIONS OF SDEs

The first theorem implies that the solution X(t) is bounded and the second theorem states that the solution is continuous on [0, T ] in the k · k^RV norm.

THEOREM: BOUNDEDNESS OF SOLUTIONS

Assume that f and g satisfy (c6) and (c7) and X ∈ H^SP is the solution of

Letting a(t) = E|X²(t)| and b(t) = 3E|X²(0)| + (3t² + 3t)k, then a(t) ≤ b(t) + (3T + 3)k

Z t 0

a(s) ds.

Using the Bellman-Gronwall inequality, it follows that a(t) ≤ b(t) + (3T + 3)k

Z t 0

exp k(3T + 3)(t − s)
b(s) ds.

As b(t) is increasing on [0, T ],

a(t) ≤ b(t) + b(t)(3T + 3)k Z t

0

exp k(3T + 3)(t − s)
ds.

Thus,

E|X(t)|² ≤ 3 E|X(0)|² + kT² + kT
exp 3k(T + T²)

for 0 ≤ t ≤ T.

THEOREM: CONTINUITY OF SOLUTIONS

Assume that f and g satisfy conditions (c6) and (c7) and X ∈ H^SP is the solution of the SDE Then,

E|X(t) − X(r)|² ≤ c|t − r| for 0 ≤ r, t ≤ T.

The previous theorem implies that there is an M > 0 such that E|X(s)|² ≤ M for 0 ≤ s ≤ T . Using this and (c7),

THE BELLMAN-GRONWALL INEQUALITY

A useful inequality is the Bellman-Gronwall inequality:

If a(t) ≤ b(t) + c

To see this, suppose that

a(t) ≤ b(t) + c

Substituting the last inequality into (*) gives:

a(t) ≤ b(t) + c Z t

0

exp(c(t − s))b(s) ds.

IT ˆO’s FORMULA AND EXACT SOLUTIONS

Itˆo’s formula is stated for stochastic differential equations.

Itˆo’s formula is helpful in finding exact solutions to certain SDEs.

Itˆo’s formula can be used to determine exact moments of the solution for certain problems.

IT ˆO’s FORMULA

Consider the Itˆo SDE in differential form:

dX(t) = f (t, X(t)) dt + g(t, X(t)) dW (t) for 0 ≤ t ≤ T with X(0) ∈ H^RV.

Let F be a smooth function so Itˆo’s formula can be applied to F (t, X). Then, F satisfies:

dF (t, X(t)) =  ∂F (t, X)

∂t + f (t, X)∂F (t, X)

∂x + 1

2g²(t, X)∂²F (t, X)

∂x²

 dt + g(t, X)∂F (t, X)

∂x dW (t).

Itˆo’s formula allows us, for example, to determine moments of the solution for certain SDEs. To find these moments, useful is the relation:

E

Z t 0

G(t, X(t)) dW (t)



= 0.

IT ˆO’s FORMULA

Before considering several examples, care must be taken with SDEs.

For example, applying Itˆo’s formula to F (t, X(t)) = X²(t) yields

d(X²(t)) = [2X(t)f (t, X(t)) + g²(t, X(t))] dt + [2X(t)g(t, X(t))] dW (t).

In particular, notice that

d(X²(t)) 6= 2X(t)dX(t) = 2X(t)[f(t, X(t)) dt + g(t, X(t)) dW (t)].

EXAMPLE: FINDING EXACT MOMENTS FOR AN SDE Consider the SDE:

dX(t) = dt + X(t) dW (t), X(0) = 0.

Then

X(t) = t + Z t

0

X(s) dW (s).

It follows that

E(X(t)) = t.

Applying Itˆo’s formula to F (t, X) = X² yields

d(X²(t)) = [2X(t) + X²(t)] dt + 2X²(t) dW (t) so that

E(X²(t)) = E Z t

0

2X(s) + X²(s)
ds.

This leads to a differential equation for E(X²(t)):

dE(X²(t))

dt = 2E(X(t)) + E(X²(t)) = 2t + E(X²(t))

with E(X²(0)) = 0. Solving this ODE gives the second moment for X(t), E(X²(t)) = −2t − 2 + 2e^t.

From these relations, Var(X(t)) = E(X²(t)) − (E(X(t))² = 2e^t − 2 − 2t − t².

This procedure can be continued. If Itˆo’s formula is applied to F (t, X) = X³, then

d(X³(t)) = [3X²(t) + 3X³(t)] dt + 3X³(t) dW (t).

This leads to the differential equation for E(X³(t)):

dE(X³(t))

dt = 3E(X²(t)) + 3E(X³(t)) = −6t − 6 + 6e^t + 3E(X³(t)) with E(X³(0)) = 0. Solving this ODE gives the third moment for X(t),

E(X³(t)) = 2t + 8

3 − 3e^t + 1 3e^3t.

EXAMPLE: FINDING EXACT MOMENTS FOR AN SDE Consider the SDE:

dX(t) = −1

4X³(t) dt + 1

2X²(t) dW (t) with X(0) = 1 2.

In this example, E(X(t)) and E(X³(t)) are to be determined exactly.

First,

dE(X(t)) = −1

4E(X³(t)) dt with E(X(0)) = 1 2

so E(X³(t)) is needed in order to find E(X(t)). Applying Itˆo’s formula to the SDE gives

dX³(t) =



−3

4X⁵(t) + 3

4X⁵(t)



dt + 3

2X⁴(t) dW (t) = 3

2X⁴(t) dW (t) with E(X³(0)) = 1

8. Thus, E(X³(t)) = 1

8 and it follows that E(X(t)) = 1

2 − 1 32t.

EXAMPLE: FINDING EXACT MOMENTS FOR AN SDE Consider the stochastic differential equation

dX(t) =  1

3X^1/3(t) + 6X^2/3(t)



dt + X^2/3(t) dW (t) with X(0) = 1.

In this example, we wish to determine E(X(t)) and E(X²(t)) exactly. First notice that

dE(X(t)) 6=  1

3 E(X(t))
1/3

+ 6 E(X(t))
2/3 dt

so an appropriate change of variables is required to find the moments. Let Y_n(t) = (X(t))^n/3 for n = 0, 1, 2, . . . , 6.

Next, applying Itˆo’s formula, the stochastic differentials are obtained dYn(t) =  1

Letting Z_n(t) = E(Y_n(t)) = E((X(t))^n/3), then the initial-value system is:

dZn(t)

dt = 1

18(n² − n)Zn−2(t) + 2nZ_n−1(t) for n = 1, 2, . . . , 6 with Zn(0) = 1 for n = 1, 2, . . . , 6 and Z0(t) = 1. Solving this gives

Z₁(t) = E((X(t))^1/3) = 2t + 1 Z₂(t) = E((X(t))^2/3) = 4t² + 37

9 t + 1 Z₃(t) = E((X(t)) = 8t³ + 38

3 t² + 19

3 t + 1 Z6(t) = E((X(t))²) = 64t⁶ + 656

3 t⁵ + 2660

9 t⁴ + 49145

243 t³ + 665

9 t² + 41

3 t + 1.

In particular, E(X(1)) = 28.0 and E(X²(1)) = 869.0206.

EXAMPLE: FINDING THE EXACT SOLUTION OF AN SDE Consider the stochastic differential equation

dX(t) = −αX(t) dt + σ dW (t), X(0) = X⁰

where α, σ, and X₀ are constants. Let F (t, X) = e^αtX(t). By Itˆo’s formula, d e^αtX(t)
= αe^αtX(t) − αe^αtX(t) dt + σe^αtdW (t)

Thus,

e^αtX(t) − X(0) = Z t

0

e^αsσ dW (s).

So the exact solution is

X(t) = X(0)e^−αt + e^−αt Z t

0

e^αsσ dW (s).

EXAMPLE: FINDING THE EXACT SOLUTION OF AN SDE Consider the stochastic differential equation

dX(t) = f (t)X(t) dt + g(t)X(t) dW (t), X(0) = X0

where X₀ is a constant. For this problem, the exact solution has the form X(t) = X₀exp

Z t

0 f (s) − 1

2g²(s)
ds + Z t

0

g(s) dW (s)

 . To see this, let F (t, X) = ln(X(t)). Applying Itˆo’s formula,

d(ln(X(t))) = f(t) − 1

2g²(t) dt + g(t) dW (t).

Thus,

ln(X(t)) − ln(X⁰) = Z t

0 f (s) − 1

2g²(s)
ds + Z t

0

g(s) dW (s) which yields the solution.

APPROXIMATING SDEs

The exact solution to an SDE is generally difficult to obtain so it is useful to be able to approximate the solution. Euler’s (or the Euler-Maruyama) method is a simple numerical method.

Euler’s method has the form

X_i+1(ω) = X_i(ω) + f (t_i, X_i(ω))∆t + g(t_i, X_i(ω))∆W_i(ω), X₀(ω) = X(0, ω)

for i = 0, 1, 2, . . . , N − 1 where Xⁱ(ω) ≈ X(tⁱ, ω), t_i = i∆t, ∆t = T /N, ∆W_i(ω) = (W (t_i+1, ω) − W (tⁱ, ω)) ∼ N(0, ∆t), and where ω indicates a sample path.

To study the error in this method, it is useful to approximate the solution for all t ∈ [0, T ] and not just at the nodal points tⁱ. To accomplish this, X(t) ≈ X(t) is defined asˆ

X(t) = Xˆ _i + Z t

t_i

f (t_i, X_i) ds + Z t

t_i

g(t_i, X_i) dW (s)

for t_i ≤ t ≤ tⁱ⁺¹ and i = 0, 1, . . . , N − 1. Notice that ˆX is identical to Euler’s method approximation at the nodal points, that is, ˆX(t_i) = X_i for i = 0, 1, . . . , N .

Notice that on the ith subinterval, ˆX(t) is the solution of the SDE

 d ˆX(t) = f (t_i, X_i) dt + g(t_i, X_i) dW (t), t_i ≤ t ≤ tⁱ⁺¹ X(tˆ _i) = X_i.

Recall that the solution X(t) satisfies the SDE

dX(t) = f (t, X(t)) dt + g(t, X(t)) dW (t), t_i ≤ t ≤ tⁱ⁺¹.

Using the inequality |2ab| ≤ a² + b² and properties of stochastic integrals, E(ǫ²(ti+1)) ≤ E(ǫ²(ti)) +

Z t_i+1

t_i E(X(t) − ˆX(t))²dt +

Z t_i+1

t_i E(f (t, X(t)) − f(tⁱ, X_i))²dt +

Z t_i+1

t_i E(g(t, X(t)) − g(tⁱ, X_i))²dt.

But

|f(t, X(t)) − f(tⁱ, X_i)|² ≤ 2|f(t, X(t)) − f(tⁱ, X(t_i))|² +2|f(tⁱ, X(t_i)) − f(tⁱ, X_i)|²

≤ 2k|t − tⁱ| + 2k|X(t) − X(tⁱ)|² + 2k|X(tⁱ) − Xⁱ|² and similarly for g using property (c6). Hence,

E(ǫ²(t_i+1)) ≤ E(ǫ²(t_i)) +

Z t_i+1

t_i E(X(t) − ˆX(t))²dt +4k(1 + c)

Z t_i+1 t_i

(t − tⁱ) dt + 4k

Z t_i+1 t_i

E(ǫ²(t_i)) dt using E|X(t) − X(tⁱ)|² ≤ c|t − tⁱ|.

Therefore,

E(ǫ²(t_i+1)) ≤ E(ǫ²(t_i))(1 + 4k∆t) + 2k(1 + c)(∆t)² +

Z t_i+1 t_i

E(ǫ²(s)) ds.

By Bellman-Gronwall inequality with b(t) = E(ǫ²(t_i))(1 + 4k∆t) +2k(1 + c)(∆t)², E(ǫ²(ti+1)) ≤ E(ǫ²(ti))(1 + 4k∆t) + 2k(1 + c)(∆t)²

+

Z t_i+1 t_i

e^(ti+1−t)E(ǫ²(t_i))(1 + 4k∆t) + 2k(1 + c)(∆t)² dt

= e^∆tE(ǫ²(t_i))(1 + 4k∆t) + 2k(1 + c)(∆t)² .

Letting ai = E(ǫ²(ti)), R = e^∆t(1 + 4k∆t), and S = e^∆t2k(1 + c)(∆t)², then a_i+1 ≤ Raⁱ + S for i = 0, 1, 2, . . . , N − 1.

These inequalities yield

a_N ≤ SR^N − 1

R − 1 with a₀ = E(ǫ²(0)) = 0.

Hence,

E(ǫ²(tN)) ≤ e^∆t2k(1 + c)(∆t)²e^{N ∆t}e^{4kN ∆t}

e^∆t − 1 + e^∆t4k∆t ≤ ∆t(1 + c)e^(1+4k)T

2 .

This holds for any nodal point and the mean square error satisfies E|X(tⁱ) − Xⁱ|² ≤ ˆc∆t

for i = 0, 1, 2, . . . , N where ˆc = ¹₂(1 + c)e^(1+4k)T.

MILSTEIN’S METHOD

Higher order numerical methods are possible for SDEs which are similar in some respects to higher order methods for ODEs. For example, there are explicit or implicit multistep methods and Runge-Kutta methods.

A popular second-order method is Milstein’s method and has mean square error proportional to (∆t)² rather than ∆t. Milstein’s method has the form

X_i+1(ω) = X_i(ω) + f (t_i, X_i(ω))∆t + g(t_i, X_i(ω))∆W_i(ω) +1

2g(t_i, X_i(ω))∂g(t_i, X_i(ω))

∂x [(∆W_i(ω))² − ∆t]

for i = 0, 1, 2, . . . , N − 1 with X⁰(ω) = X(0, ω), where X_i(ω) ≈ X(tⁱ, ω), ∆W_i(ω) = (W (ti+1, ω) − W (tⁱ, ω)) ∼ N(0, ∆t), tⁱ = i∆t, ∆t = T /N, and where ω indicates a sample path.

Milstein’s method has an additional term at each step in comparison with Euler’s method.

EXAMPLE: APPROXIMATION OF AN SDE BY EULER AND MIL-STEIN

Consider the stochastic differential equation dX(t) =  1

3X^1/3(t) + 6X^2/3(t)



dt + X^2/3(t) dW (t), X(0) = 1.

It was shown earlier that E(X(1)) = 28.0 and E(X²(1)) = 869.0206. For this problem, Euler’s method has the form:

X_i+1 = X_i + 1

3X_i^1/3 + 6X_i^2/3



∆t + X_i^2/3√

∆t η_i where η_i ∼ N(0, 1)

for i = 0, 1, 2, . . . , N − 1 with X⁰ = 1, ti = i∆t, and ∆t = 1/N . Milstein’s method has the form

X_i+1 = X_i + 1

3X_i^1/3 + 6X_i^2/3



∆t + X_i^2/3√

∆t η_i + 1

3X_i^1/3(η_i² − 1)∆t where ηi ∼ N(0, 1).

The calculational results for the mean square error E|X(1) − X^N|² are given in the table for 10,000 sample paths for each value of N .

Value of N Euler Error Milstein Error 2⁹ 2.80× 10⁻² 1.61× 10⁻²

2¹⁰ 1.04× 10⁻² 4.03× 10⁻³ 2¹¹ 4.20× 10⁻³ 1.01× 10⁻³ 2¹² 1.89× 10⁻³ 2.53× 10⁻⁴ 2¹³ 8.76× 10⁻⁴ 6.24× 10⁻⁵ 2¹⁴ 4.12× 10⁻⁴ 1.60× 10⁻⁵

Notice that the mean square errors are approximately proportional to ∆t = 1/N for Euler’s method and to (∆t)² = 1/N² for Milstein’s method.

Next, for this example, E(X(1)) and E(X(1))² were estimated using E(X(1)) ≈ P100,000

j=1 X_N^(j)/100, 000 and E(X(1))² ≈ P100,000

j=1 (X_N^(j))²/100, 000 where X_N^(j) is the es-timate of X(1) for the jth sample path using N intervals.

In the table, the errors are given in parentheses. Recall that E(X(1)) = 28.0 and E(X²(1)) = 869.0206 are the exact values. Notice that the errors in the mean values are proportional to ∆t for either numerical method. In particular, the errors in Euler’s method when estimating mean values are proportional to ∆t rather than (∆t)^1/2.

Value of N Euler Estimate Milstein Estimate 2⁶ 27.07 (0.93) 27.08 (0.92)

2⁷ 27.56 (0.44) 27.56 (0.44) 2⁸ 27.79 (0.21) 27.79 (0.21)

Value of N Euler Estimate Milstein Estimate 2⁶ 810.15 (58.87) 810.18 (58.84) 2⁷ 840.89 (28.13) 840.93 (28.09) 2⁸ 855.33 (13.69) 855.31 (13.71)

0 0.2 0.4 0.6 0.8 1 0

5 10 15 20 25 30

Time

Mean and One Sample Path

Figure 1: Mean solution and one sample path

The mean and one sample path are plotted in the figure for this problem.

STRONG AND WEAK APPROXIMATIONS

This example illustrates that are two kinds of approximation commonly discussed in computational solution of SDEs.

A numerical method is said to be a strong approximation of order γ if kX(T ) − X^Nk^RV ≤ c(∆t)^γ

where X(T ) is the exact solution at time T and X_N is the approximate solution using step length ∆t = T /N .

Euler’s method and Milstein’s methods have strong orders ¹₂ and 1.

However, if expectations of functions of a solution to an SDE are desired and not necessarily the pathwise approximation provided by a strong ap-proximation, then a weak numerical method may be sufficient.

An approximation X_N is said to converge weakly with order β if

|E(F (X(T ))) − E(F (X^N))| ≤ c(∆t)^β

for all smooth functions F , where ∆t = T /N is the step size.

Euler’s method and Milstein’s method both have weak order 1.

RICHARDSON EXTRAPOLATION

Both Euler’s or Milstein’s method have weak-error expansions of the cor-rect form for applying Richardson extrapolation.

The weak error for Euler’s or Milstein’s method has the form E(F (X(T ))) − E(F (X^N)) = c₁∆t + c₂(∆t)² + c₃(∆t)³ + . . . , where ∆t = T /N and c₁, c₂, c₃, . . . are independent of ∆t.

So, several approximations with different values of N can be applied to obtain a higher order approximation. Suppose that E(F (XN)), E(F (X2N)), and E(F (X4N)) are three approximations to E(F (X(T ))) using step lengths of T /N , T /2N , and T /4N in Euler’s method or in Milstein’s method.

To obtain an approximation to E(F (X(T ))) of order (∆t)², let

E(F (X(T ))) − 2E(F (X_2N)) − E(F (X^N)) = ˆc₂(∆t)² + ˆc₃(∆t)³ + . . . . To obtain an approximation to E(F (X(T ))) of order (∆t)³, let

E(F (X(T ))) −8E(F (X4N)) − 6E(F (X^2N)) + E(F (XN))/3 = ˜c3(∆t)³ + . . . .

EXAMPLE: RICHARDSON EXTRAPOLATION

Referring to the values in for the previous example, the following approxi-mations to E((X(1))²) are obtained using Euler’s method:

E((X₆₄)²) = 810.15, E((X₁₂₈)²) = 840.89, and E((X₂₅₆)²) = 855.33.

To obtain O((∆t)²) and O((∆t)³) approximations, respectively, to E((X(1))²) we calculate

2E((X128)²) − E((X⁶⁴)²) = 871.63 and

[8E((X₂₅₆)²) − 6E((X¹²⁸)²) + E((X₆₄)²)]/3 = 869.15.

As E((X(1))²) = 869.02 exactly, the original Euler approximations are much improved through extrapolation.

STRONG APPROXIMATIONS ARE WEAK APPROXIMATIONS

Any strong approximation is also a weak approximation as the Lyapunov inequality gives

|E(F (X(T ))) − E(F (X^N))| ≤ (E|(F (X(T ))) − (F (X^N))|²)^1/2

≤ L(E|X(T ) − X^N|²)^1/2 = LkX(t) − X^Nk^RV assuming that F satisfies a Lipschitz condition.

However, there are weak methods which are not strong approximations.

Consider the discrete process described earlier. For a particular trajectory, suppose that at time t_i, X_i = mδ for some integer m, where δ > 0 is small.

Define the three possibilities at time t_i+1 = t_i + ∆t as







X_i+1 = X_i + δ with probability r(t_i, X_i)∆t/δ²,

Xi+1 = Xi with probability 1 − r(tⁱ, Xi)∆t/δ² − s(tⁱ, Xi)∆t/δ², Xi+1 = Xi − δ with probability s(tⁱ, Xi)∆t/δ²,

where

 r(t_i, X_i) = f (t_i, X_i)δ + g²(t_i, X_i)
/2 s(t_i, X_i) = − f(tⁱ, X_i)δ + g²(t_i, X_i)
/2.

The probability distribution of XN approaches that of X(T ) as ∆t, δ → 0 implying that E(F (X_N)) ≈ E(F (X(T ))) for small values of ∆t and δ.

Apply the weak method to the previous problem. The calculational re-sults are presented using 100,000 sample paths for estimating E(X(1)) and E(X(1))². Recalling that E(X(1)) = 28.00 and E(X(1))² = 869.02, the calcula-tional results are reasonable.

Value of N E(X(1)) Estimate E(X(1))² Estimate

2¹² 17.04 292.42

2¹³ 24.64 630.81

2¹⁴ 27.88 854.12

2¹⁵ 27.97 868.12

SYSTEMS OF STOCHASTIC DIFFERENTIAL EQUATIONS

Itˆo’s formula and numerical methods can be extended to systems. Let X(t, ω) = [X~ 1(t, ω), X2(t, ω), . . . , Xd(t, ω)]^T

W (t, ω) = [W~ ₁(t, ω), W₂(t, ω), . . . , W_m(t, ω)]^T f : [0, T ] × R~ ^d → R^d and g : [0, T ] × R^d → R^d×m, where Wi(t, ω), 1 ≤ i ≤ m are independent Wiener processes.

Then a system of stochastic differential equations has the form d ~X(t, ω) = ~f (t, ~X(t, ω)) dt + g(t, ~X(t, ω)) d ~W (t, ω).

In component form, the system is X_i(t) = X_i(0) +

Z t 0

f_i(s, ~X(s)) ds +

m

X

j=1

Z t 0

g_i,j(s, ~X(s)) dW_j(s) for i = 1, 2, . . . , d.

IT ˆO’s FORMULA FOR SYSTEMS Let

F : [0, T ] × R~ ^d → R^k and let ~Y (t, ω) = ~F (t, ~X(t, ω)).

Then the pth component of ~Y (t, ω) satisfies:

dY_p(t) =

EXAMPLE: IT ˆO’s FORMULA FOR A PROBLEM WITH d = 1 AND m = 2 Consider the SDE:

 dX(t) = t²X(t) dt + t dW1(t) + X(t) dW2(t), 0 ≤ t ≤ T X(0) = 1,

where d = 1 and m = 2.

For this problem, f₁ = t²X, g_1,1 = t, and g_1,2 = X.

Consider using Itˆo’s formula to find the SDE for F = X². Applying Itˆo’s formula,

 d(X²(t)) = 2t²X²(t) + t² + X²(t) dt + 2tX(t) dW1(t) + 2X²(t) dW₂(t) X²(0) = 1.

EULER’S AND MILSTEIN’S METHODS FOR SYSTEMS Euler’s method for systems has the form

X~n+1(ω) = ~Xn(ω) + ~f (tn, ~Xn(ω))∆t + g(tn, ~Xn(ω))∆ ~Wn(ω)

for n = 0, 1, 2, . . . , N , where ~X_n(ω) ≈ ~X(t_n, ω), ∆t = T /N , ∆ ~W_n = ~W (t_n+1)− ~W (t_n).

In component form, Euler’s method is:

X_i,n+1(ω) = X_i,n(ω) + f_i(t_n, ~X_n(ω))∆t +

m

X

j=1

g_i,j(t_n, ~X_n(ω))∆W_j,n(ω) for i = 1, 2, . . . , d, where ∆W_j,n ∼ N(0, ∆t).

Milstein’s method for multidimensional SDEs involves the double stochastic integral

I_n(j₁, j₂) =

Z tn+∆t t_n

Z s t_n

dW_j1(r) dW_j2(s).

Milstein’s method has the componentwise form X_i,n+1(ω) = X_i,n(ω) + f_i(t_n, ~X_n(ω))∆t +

m

X

j=1

g_i,j(t_n, ~X_n(ω))∆W_j,n(ω)

+

m

X

j₁=1 m

X

j₂=1 d

X

l=1

g_l,j1∂gi,j₂

∂x_l I_n(j₁, j₂) for i = 1, 2, . . . , d.

EXAMPLE: APPROXIMATION OF AN SDE WITH d = 1 AND m = 2 Consider the SDE;

 dX(t) = t²X(t)dt + tdW₁(t) + X(t)dW₂(t), 0 ≤ t ≤ T X(0) = 1,

where d = 1 and m = 2.

For this problem, Euler’s method has the form

 Xn+1 = Xn + t²_nXn∆t + tn∆W1,n + Xn∆W2,n

X₀ = 1,

for n = 0, 1, 2, . . . , where ∆W_1,n, ∆W_2,n ∼ N(0, ∆t) and tⁿ = n∆t.

Milstein’s method has the form

 X_n+1 = X_n + t²_nX_n∆t + t_n∆W_1,n + X_n∆W_2,n + t_nI_n(1, 2) + X_nI_n(2, 2) X₀ = 1,

for n = 0, 1, 2, . . . .

It is useful to note that I_n(j₁, j₁) =

Z tn+∆t tn

Z s tn

dW_j1(r) dW_j1(s) = 1

2 (∆W_j1,n)² − ∆t

but I_n(j₁, j₂) for j₁ 6= j² does not have an analytical form and must be ap-proximated.

This multiple integral can be approximated by a Fourier series expansion.

Also, if [tn, tn+1] is divided into M equal intervals with tj,n = tn + j∆t/M for j = 0, 1, . . . , M, then

I_n(j₁, j₂) ≈ ˜I_n(j₁, j₂) =

M −1X

j=0

[W_j1(t_j,n)− W^j1(t_0,n)][W_j2(t_j+1,n) − W^j2(t_j,n)].

It can be shown that E|Iⁿ − ˜I_n|² = (∆t)²/(2M ).

FORWARD KOLMOGOROV (FOKKER-PLANCK) EQUATION

The probability distribution of solutions to a discrete-valued continuous stochastic process satisfies a system of differential equations called the for-ward Kolmogorov equations. An analogous result holds for the probability distribution of solutions to an SDE.

Consider the stochastic differential equation

dX(t) = f (t, X(t)) dt + g(t, X(t)) dW (t)

Let p(t, x) be the probability density for solutions to the SDE. The previous

Integrating by parts the right-hand side yields the relation Z _∞

As the above integral holds for every function F ∈ C^∞0 (R), this implies that

∂p(t, x)

This equation is the forward Kolmogorov equation or Fokker-Planck equa-tion for the probability distribuequa-tion of soluequa-tions to the SDE.

Also, the forward Kolmogorov equation for a system of SDEs is:

∂p(t, ~x)

EXAMPLE: SOLUTION OF A FORWARD KOLMOGOROV EQUATION Consider the stochastic differential equation

 dX(t) = a dt + b dW (t) X(0) = x₀.

The probability density of the solutions satisfies the forward Kolmogorov equation







∂p(t, x)

∂t = −∂(ap(t, x))

∂x + b² 2

∂²(p(t, x))

∂²x p(0, x) = δ(x − x⁰).

The solution to this partial differential equation is p(t, x) = 1

√2πb²t exp −(x − at − x⁰)² 2b²t

 .

STABILITY

There are several kinds of stability questions and several ways to define stability for SDEs. To introduce this topic, it is useful to first review stability concepts for ODEs. Consider the initial-value problem:







d~y(t)

dt = ~f (~y(t)) for t > 0

~y(0) = ~a

where ~y : R → Rⁿ and ~f : Rⁿ → Rⁿ. Suppose that ~z(t) satisfies the same differential equation as ~y(t) but with a different initial condition, i.e.,







d~z(t)

dt = ~f (~z) for t > 0

~z(0) 6= ~a.

Suppose that ~a = ~γ is a critical point, i.e., ~f (~γ) = ~0. Then the solution satisfies ~y(t) = ~γ for t ≥ 0. The initial-value problem is said to be stable at

~γ if given ǫ > 0 there is a δ > 0 such that

k~z(t) − ~γk < ǫ for t ≥ 0 whenever k~z(0) − ~γk < δ.

That is, small changes in the initial condition do not produce large changes in the solution for t ≥ 0.

In numerical solution of initial-value problems for ODEs, there are two common numerical stability concepts. Suppose that a single-step method for solving has the form:

 ~y_k+1 = ~y_k + h ~φ(h, ~y_k) for k = 0, 1, 2, . . . , N − 1

~y₀ = ~a.

where h = T /N is the step length, tk = kh, and ~yk ≈ ~y(t^k) for each 0 ≤ k ≤ N.

The method is numerically stable if small changes in the initial condition do not produce large changes in the computational solution.

Specifically, if ~z_k for k = 0, 1, . . . , N satisfies the method but with a different initial condition ~z₀ 6= ~y⁰, then the numerical scheme is numerically stable provided that

k~y^k − ~z^kk ≤ cǫ for 0 ≤ k ≤ N when k~y⁰ − ~z⁰k < ǫ.

If ~φ satisfies an appropriate Lipschitz condition, then the numerical scheme can be shown to be stable. However, the constant c can be extremely large, especially for stiff systems, which motivates another concept of numerical stability.

To study stability of stiff systems, the following scalar test problem is

stud-ied: 





dy(t)

dt = λy for t > 0 y(0) = a

where λ is a constant. Clearly, if a 6= 0, then y(t) → 0 as t → ∞ if and only if Re(λ) < 0.

Consider, for example, applying Euler’s method to this test problem. Then

 yk+1 = (1 + hλ)yk for k = 0, 1, 2, . . . , y₀ = a.

and y_k → 0 as k → ∞, if and only if −2 < Re(λh) < 0. The region of absolute stability of Euler’s method is −2 < Re(λh) < 0. The region of absolute stability gives a condition on the step length. If the method satisfies this condition, then the numerical solution does not “blow up” but decreases to zero behaving like the solution to the initial-value problem. For Euler’s method to behave similarly to the solution for a large negative value of λ, the step length h must be selected to be small.

However, for the backward Euler method, which for the test problem has the form:

 y_k+1 = y_k + hλy_k+1 for k = 0, 1, 2, . . . , y₀ = a,

the region of absolute stability is the entire left-half of the complex plane, i.e. −∞ < Re(λh) < 0.

To see this, consider

y_k+1 = y_k/(1 − hλ) for k = 0, 1, 2, . . . .

For this method, the step length h need not be chosen very small for the numerical solution to perform similarly to the actual solution even for an initial-value problem that involves a large negative value of λ.

The concept of absolute stability is useful when considering numerical so-lution of systems. For the test initial-value system







d~y(t)

dt = A~y for 0 ≤ t ≤ T

~y(0) = ~a.

where A is an n × n matrix, then ~y(t) → ~0 as t → ∞ provided that Re(λⁱ) < 0 for each eigenvalue λ_i for 1 ≤ i ≤ n. For this problem, Euler’s method is

 ~yk+1 = (I + Ah)~yk for k = 0, 1, 2, . . . ,

~y₀ = ~a.

The eigenvalues of I + Ah are 1 + λih for i = 1, 2, . . . , n and ~yk → ~0 as k →

∞ provided that −2 < Re(λⁱ)h < 0 for each eigenvalue λi for 1 ≤ i ≤ n.

Hence, the step length h is forced to satisfy a condition determined by the eigenvalues with large negative real parts.

Now consider stability for SDEs. First, stability of a steady solution to an SDE is studied then numerical stability of an approximation is studied.

Consider stability of a steady solution for the SDE

 dX(t) = f(X(t)) dt + g(X(t)) dW (t) for 0 ≤ t ≤ T X(0) = a.

It is supposed that f (0) = g(0) = 0 so that X(t) ≡ 0 is a steady solution.

There are many ways to define stochastic stability for a steady solution of an SDE. Two ways are considered here, asymptotic stochastic stability and mean-square stability. It is assumed that X(0) 6= 0.

If lim

t→∞|X(t)| = 0 with probability 1, then X(t) ≡ 0 is said to be asymptotically stochastically stable.

If lim

t→∞E(|X(t)|²) = 0, then X(t) ≡ 0 is said to be mean-square stable.

It is interesting that some SDEs may be both asymptotically stochastically stable and mean-square stable while others may be asymptotically stochas-tically stable but not mean-square stable.

To illustrate this behavior, stability is analyzed for an SDE with f (X) = λX and g(X) = µX. In this case,

 dX(t) = λX(t) dt + µX(t) dW (t) for 0 ≤ t ≤ T X(0) = a

and E(X(t)) = X(0) exp(λt).

Using Itˆo’s formula, X²(t) satisfies the SDE

 d(X²(t)) = 2λX²(t) + µ²X²(t)
dt + 2µX²(t) dW (t) for t > 0 X²(0) = a².

It follows that E(X²(t)) satisfies the SDE

 d(E(X²(t))) = 2λE(X²(t)) + µ²E(X²(t))
dt for t > 0 E(X²(0)) = a².

and the solution E(X²(t)) is found to be

E(X²(t)) = X²(0) exp((2λ + µ²)t).

This solution implies that the steady solution X(t) = 0 is mean-square stable if and only if λ + µ²/2 < 0.

Now consider Itˆo’s formula applied to ln(X(t)). Then,

 d(ln(X(t)) = (λ − µ²/2) dt + µ dW (t) for t > 0 ln(X(0)) = ln(a).

Let ∆t be a given interval width and let t_i = i∆t for t = 0, 1, 2, . . . . This SDE can be exactly integrated from ti to ti+1 to yield:

ln(X(ti+1)) − ln(X(tⁱ)) = (λ − µ²/2) (ti+1 − tⁱ) + µηi p(t_i+1 − tⁱ)

. By the Law of Large Numbers, S_n

But, letting t = tn, 1 n∆t

Xn−1 i=0

ln X(t_i+1) X(t_i)



= 1

n∆t ln X(t_n) X(0)



= 1

t ln X(t) X(0)



→ (λ − µ²/2) w.p.1 as t → ∞.

Therefore,

X(t) → X(0) exp((λ − µ²/2)t) w.p.1 as t → ∞.

This result implies that the steady solution X(t) = 0 is asymptotically sto-chastically stable if and only if λ − µ²/2 < 0.

Hence, for example, if λ = µ²/4 in the SDE, then X(t) → 0 with probability 1 as t → ∞ while E(X(t)) → ∞ and E(X²(t)) → ∞ under the same condition.

Now, numerical stability of SDEs is considered, in particular, with respect to stiff stochastic problems with additive noise and then, more briefly, with respect to multiplicative noise.

The test problem for additive noise has the form

 dX(t) = λX(t) dt + µ dW (t) for t > 0 X(0) = a.

Two kinds of numerical stochastic stability are numerical asymptotic sto-chastic stability and numerical mean-square stability. Let X_k and ˜X_k be two approximations with the same numerical method but with different initial values.

If lim

k→∞|X^k − ˜X_k| = 0 with probability 1, then the approximation is said to be asymptotically stochastically stable.

If lim

k→∞E(|X^k − ˜Xk|²) = 0, then the approximation is said to be mean-square stable.

Consider first Euler’s method for solution of test problem:

 X_k+1 = X_k + λX_kh + µ η_k√

h for k = 0, 1, . . . X₀ = a

where Xk ≈ X(kh), η^k ∼ N(0, 1) for each k, and h is the step length. Fur-thermore, let ˜X_k be another numerical approximation but with a different initial approximation ˜X₀ = ˜a. Let Z_k = X_k − ˜X_k. Then, Z_k satisfies

 Z_k+1 = Z_k + λh Z_k for k = 0, 1, . . . Z0 = a − ˜a

and therefore,

|X^k − ˜X_k| = |Z^k| = |1 + λh|^k |Z⁰| for k = 0, 1, . . . .

Thus, Euler’s method is asymptotically and mean square stable for the test problem provided that −2 < λh < 0. An analogous result holds for stability of Euler’s method for stiff systems with additive noise. Specifically, Euler’s method is numerically stable for a stochastic system with additive noise

( d ~X(t) = A ~X(t) dt + µ d ~W (t) for t > 0 X(0) = ~a.~

provided that −2 < Re(λⁱ)h < 0 for each eigenvalue λ_i of A.

在文檔中 An Intensive Course in Modeling Techniques and Numerical Methods for Stochastic Differential Equations (頁 102-200)