Efficient Variable Step Size Approximations for Strong Solutions of Stochastic Differential Equations with Additive Noise and Time Singularity

Abstract

We consider stochastic differential equations with additive noise and conditions on the coefficients in those equations that allow a time singularity in the drift coefficient. Given a maximum step size, ℎ^∗, we specify variable (adaptive) step sizes relative to ℎ^∗ which decrease as the time node points approach the singularity. We use an Euler-type numerical scheme to produce an approximate solution and estimate the error in the approximation. When the solution is restricted to a fixed closed time interval excluding the singularity, we obtain a global pointwise error of order 𝑂(ℎ^∗). An order of error 𝑂(ℎ^∗𝑝) for any 𝑝<1 is obtained when the approximation is run up to a time within ℎ^∗𝑞 of the singularity for an appropriate choice of exponent 𝑞. We apply this scheme to Brownian bridge, which is defined as the nonanticipating solution of a stochastic differential equation of the type under consideration. In this special case, we show that the global pointwise error is of order 𝑂(ℎ^∗), independent of how close to the singularity the approximation is considered.

1. Introduction

Numerical approximation methods for stochastic differential equations (SDEs) are well developed for SDEs with coefficients satisfying Lipschitz conditions. One of basic methods is the Euler-Maruyama algorithm which is presented in [1] and also well explained by Kloeden and Platen [2]. It is developed using Taylor series expansions and the Itô formula and is analogous to the Euler algorithm for the deterministic case. The Euler method with constant step sizes gives the global error 𝑂(ℎ^1/2) when drift and diffusion terms satisfy global Lipschitz and linear growth conditions. In addition to methods with constant step sizes, methods with variable (adaptive) step sizes have been developed, for example, those in [3–7]. In general, these methods seek to improve the efficiency of the algorithms by adapting the step sizes to the state of the numerical solution at successive stages of the algorithm.

We consider stochastic differential equations with additive noise and an endpoint singularity with respect to the time variable in the drift term. With this time singularity, the global error estimate for the Euler-Maruyama scheme can only be applied up to a fixed time before the singularity. Using fixed step sizes to generate a numerical approximation for a time closer to the singularity comes at a great cost in efficiency due to the much smaller step size required to produce the same order of global error. We seek to increase the efficiency of the algorithm by using variable step sizes adapted to the shape of the singularity. In making this adaptation, we must also ensure that the step sizes do not become so small that the sum of the steps cannot reach times arbitrarily close to the singularity.

We consider a stochastic differential equation of the form𝑑𝑋_𝑡=𝑔𝑡,𝑋_𝑡1−𝑡𝑑𝑡+𝜎𝑑𝑊_𝑡,0≤𝑡<1,𝑋(0)=𝑥₀, (1a)(1b)where 𝜎 is a constant, 𝑊_𝑡 is standard Brownian motion, and 𝑔(𝑡,𝑥) satisfies linear growth and Lipschitz conditions with bounds for derivatives up to the second order. For example, suppose that, uniformly in 𝑡 on [0,1),𝑔(𝑡,𝑥)≤𝑘₀(1+|𝑥|),𝑔₁(𝑡,𝑥)≤𝑘₁(1+|𝑥|),𝑔₂(𝑡,𝑥)≤𝐶₁,𝑔₂₂(𝑡,𝑥)≤𝐶₂,𝔼𝑋²_𝑡≤𝐾, (2a)(2b)(2c)(2d)(2e)where 𝑔_𝑖 denotes the partial derivative of 𝑔 with respect to the 𝑖th variable.

Note that the singularity in the drift term stops us from using the standard results of numerical analysis. To overcome this problem, we use variable step sizes and stop the approximation at 𝑡=1−𝛿 where 0<𝛿<1. We estimate bounds for the pointwise local error and we use a lemma [8] to estimate the 𝐿² global error. We fix a maximum step size ℎ^∗ and produce variable steps and node points depending on ℎ^∗ in order to produce an estimate of the global error of order 𝑂(ℎ^∗) for fixed 𝛿. The same algorithm and analysis using fixed step sizes fail to produce an estimate of global error. We further show that there is a choice of 𝛿 depending on ℎ^∗ so that the global error is of order 𝑂(ℎ^∗𝑝) for positive 𝑝 and 𝛿 converges to 0 as ℎ^∗ converges to 0. In this analysis we also obtain an estimate of the number of node points required for a given ℎ^∗ and 𝛿. The trade-off for producing our global error estimate is that this number grows faster than 1/ℎ^∗, the corresponding number of node points for fixed step size approximations.

In Section 4, we apply this method to the process well known as the Brownian bridge (BB). Here we use the version of Brownian bridge that is the strong, adapted solution of a stochastic differential equation of the form (1a)-(1b) [9, p. 75]. BB is a stochastic process that starts at a specific point and converges to a specific point at a given time. BB is widely used in Monte Carlo simulations and many other applications [10, 11]. Most of the applications use simpler anticipating versions of BB like 𝑋_𝑡=𝑊_𝑡−𝑡𝑊₁ to avoid the approximation problem that occurs from the singularity in the drift term. We apply our methods to numerically approximate Brownian bridge and obtain even better estimates of global error than those for our general problem. For more properties of Brownian bridge, one can refer to [12, p. 86–89].

Brownian bridge from 𝑎 to 𝑏 is given by the following SDE [9, p. 75]. For the sake of simplicity, we use 𝑎=0 and 𝑏=1 later in our analysis. For fixed 𝑎, 𝑏∈ℝ,𝑑𝑋_𝑡=𝑏−𝑋_𝑡1−𝑡𝑑𝑡+𝑑𝑊_𝑡,0≤𝑡<1,𝑋₀=𝑎. (3a)(3b)Then the solution of the above SDE is given by 𝑋_𝑡=𝑎(1−𝑡)+𝑏𝑡+(1−𝑡)^𝑡₀𝑑𝑊_𝑠1−𝑠,0≤𝑡<1. (4) Also, lim_𝑡→1𝑋_𝑡=𝑏,almostsurely. (5) It is clear that BB satisfies conditions given in (2a)–(2d). It can also be proven that BB satisfies the bounded moment condition, (2e), but we will show that our proposed numerical method can be applied directly to the Brownian bridge process to obtain the improved global error estimate without using the assumption of bounded moments.

In Section 2, we discuss the Euler-Maruyama method and some lemmas used in this paper. In Section 3, we propose the numerical approximation with variable (adaptive) step sizes based on the Euler-Maruyama algorithm and we analyze the error. In Section 4, we apply these methods to approximate the specific example, Brownian bridge, and we show that the order of the global error estimate is better for BB than for the general case.

2. Euler-Maruyama Algorithm and Fundamental Lemmas

Let 𝑊_𝑡 be standard Brownian motion, and let ℱ_𝑡 be the natural filtration. We use the following general stochastic differential equation and Itô formula to present the development of Euler-Maruyama algorithm: 𝑑𝑋_𝑡=𝑎𝑡,𝑋_𝑡𝑑𝑡+𝑏𝑡,𝑋_𝑡𝑑𝑊_𝑡,𝑋₀=𝜉, (6) where 𝑎(𝑡,𝑋_𝑡),𝑏(𝑡,𝑋_𝑡)∈[0,∞]×Ω→ℝ.

The same thing can be given in integral notation assuming all the integrals exist [9, p. 22]. The second integral is the Itô integral 𝑋_𝑡=𝑋₀+^𝑡₀𝑎𝑠,𝑋_𝑠𝑑𝑠+^𝑡₀𝑏𝑠,𝑋_𝑠𝑑𝑊_𝑠. (7)

Theorem 1 (Itô Formula [9, p. 46]). Let 𝑋_𝑡 be an Itô process given by 𝑑𝑋_𝑡=𝑎𝑑𝑡+𝑏𝑑𝑊_𝑡, (8) where 𝑎∶[0,∞]×Ω→ℝ and 𝑏∶[0,∞]×Ω→ℝ are measurable, ℱ_𝑡-adapted, and integrable (resp., square integrable) on every interval [0,𝑡], almost surely. Let 𝑔(𝑡,𝑥)∈𝒞²([0,∞)×ℝ). Then 𝑌_𝑡=𝑔(𝑡,𝑋_𝑡) is again an Itô process and 𝑑𝑌_𝑡=𝜕𝑔𝜕𝑡𝑡,𝑋_𝑡𝑑𝑡+𝜕𝑔𝜕𝑥𝑡,𝑋_𝑡𝑑𝑋_𝑡+12𝑏²𝜕²𝑔𝜕𝑥²𝑡,𝑋_𝑡𝑑𝑡. (9)

Consider a sequence of times 𝑡₀<𝑡₁<𝑡₂<⋯<𝑡_𝑁 and the corresponding step sizes ℎ_𝑛=𝑡_𝑛−𝑡_𝑛−1. Following [13, p. 56], we use the Itô formula to expand the drift and diffusion terms in (7), 𝑋𝑡_𝑛+1=𝑋𝑡_𝑛+^𝑡𝑛+1_𝑡𝑛𝑎_𝑡𝑑𝑡+^𝑡𝑛+1_𝑡𝑛𝑏_𝑡𝑑𝑊_𝑡=𝑋𝑡_𝑛+^𝑡𝑛+1_𝑡𝑛𝑎_𝑡𝑛+^𝑡_𝑡𝑛(𝐴𝑎)_𝑠𝑑𝑠+^𝑡_𝑡𝑛(𝐵𝑎)_𝑠𝑑𝑊_𝑠+𝑑𝑡_{𝑡𝑛+1𝑡𝑛}𝑏_𝑡𝑛+^𝑡_𝑡𝑛(𝐴𝑏)_𝑠𝑑𝑠+^𝑡_𝑡𝑛(𝐵𝑏)_𝑠𝑑𝑊_𝑠𝑑𝑊_𝑡=𝑋𝑡_𝑛+𝑎_𝑡𝑛ℎ_𝑛+1+𝑏_𝑡𝑛Δ𝑊_𝑛+1+^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛(𝐴𝑎)_𝑠𝑑𝑠𝑑𝑡+^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛(𝐵𝑎)_𝑠𝑑𝑊_𝑠𝑑𝑡+^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛(𝐴𝑏)_𝑠𝑑𝑠𝑑𝑊_𝑡+^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛(𝐵𝑏)_𝑠𝑑𝑊_𝑠𝑑𝑊_𝑡, (10) where ℎ_𝑛+1=𝑡_𝑛+1−𝑡_𝑛, Δ𝑊_𝑛+1=𝑊(𝑡_𝑛+1)−𝑊(𝑡_𝑛), 𝐴𝑓=𝑓₁+𝑎𝑓₂+(1/2)𝑏²𝑓₂₂, 𝐵𝑓=𝑏𝑓₂, 𝑎_𝑠=𝑎(𝑠,𝑋(𝑠)), and 𝑏_𝑠=𝑏(𝑠,𝑋(𝑠)). Because 𝑊 is standard Brownian motion, Δ𝑊_𝑛 forms a sequence of independent Gaussian random variables with mean zero and Var(Δ𝑊_𝑛)=ℎ_𝑛.

The first three terms of expansion generate the Euler-Maruyama algorithm with variable step sizes ℎ_𝑛 [2, p. 190]. The double integrals yield the local truncation error. So we define the numerical scheme by 𝑌_𝑛+1=𝑌_𝑛+𝑎𝑡_𝑛,𝑌_𝑛ℎ_𝑛+1+𝑏𝑡_𝑛,𝑌_𝑛Δ𝑊_𝑛+1. (11)

With suitable conditions on the coefficient function 𝑎 and with 𝑏 constant, it can be shown that mean square local truncation error is of order 𝑂(ℎ³_𝑛+1) as ℎ_𝑛+1↓0. Let 𝜖_𝑛+1=𝑋_𝑛+1−𝑌_𝑛+1. When step sizes are fixed at ℎ, it follows that [13, p. 60] 𝔼𝜖²_𝑛+1≤1+𝑚₁ℎ𝔼𝜖²_𝑛+𝑚₂ℎ³. (12) It can be further shown that the global error (𝔼[𝜖²_𝑛])^1/2=𝑂(ℎ) as ℎ↓0 [13, p. 61]. In the process of estimating the global error, the following lemma is used. This is an easy extension of a fundamental lemma of numerical analysis [8, p. 189].

Lemma 2. If 𝑀₁≥−1 and 𝑀₂≥0 are real numbers and {𝑎_𝑛}^𝑁_𝑛=0 is a sequence with 𝑎₀≥0 s.t. 𝑎_𝑛+1≤1+𝑀₁𝑎_𝑛+𝑀₂ (13) for all 𝑛=0,1,2,…,𝑁−1, then 𝑎_𝑛+1≤𝑒^𝑛𝑀1𝑀₂𝑀₁+𝑎₀−𝑀₂𝑀₁ (14) for all 𝑛=0,1,2,…,𝑁−1.

3. Numerical Approximation for SDE with Singularity

We return to the SDE of the form (1a)-(1b) presented in Section 1:𝑑𝑋_𝑡=𝑔𝑡,𝑋_𝑡1−𝑡𝑑𝑡+𝜎𝑑𝑊_𝑡,0≤𝑡<1,𝑋(0)=𝑥₀, (15a)(15b)where 𝜎 is a constant and 𝑔(𝑡,𝑥) satisfies linear growth and Lipschitz conditions with bounds for derivatives up to the second order. In particular, we suppose that, uniformly in 𝑡 on [0,1),𝔼𝑔𝑡,𝑋_𝑡²≤𝐾₀,𝔼𝑔₁𝑡,𝑋_𝑡²≤𝐾₁,𝑔₂(𝑡,𝑥)≤𝐶₁,𝑔₂₂(𝑡,𝑥)≤𝐶₂. (16a)(16b)(16c)(16d)It follows then that 𝑔(𝑡,𝑥)−𝑔𝑡,𝑦≤𝐶₁𝑥−𝑦. (17) These conditions follow, for example, from the conditions given in (2a)–(2e).

In this case, 𝑎(𝑡,𝑥)=𝑔(𝑡,𝑥)/(1−𝑡) and 𝑏(𝑡,𝑥)=𝜎, and thus 𝑔()𝐴𝑎=𝑡,𝑥(1−𝑡)²+𝑔₁(𝑡,𝑥)+𝑔()𝑔1−𝑡𝑡,𝑥₂(𝑡,𝑥)(1−𝑡)²+𝜎²𝑔₂₂(𝑡,𝑥)2(),1−𝑡𝐵𝑎=𝜎𝑔₂(𝑡,𝑥)(),1−𝑡𝐴𝑏=𝐵𝑏=0. (18) Note that the derivatives of 𝑎 are not bounded. Nevertheless, on each subinterval where 𝑡_𝑛+1<1, we have the expansion from (10) and the numerical algorithm from (11), respectively: 𝑋_𝑛+1=𝑋_𝑛+𝑔𝑡_𝑛,𝑋_𝑛1−𝑡_𝑛ℎ_𝑛+1+𝜎Δ𝑊_𝑛+1+^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛(𝐴𝑎)_𝑠𝑑𝑠𝑑𝑡+^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛(𝐵𝑎)_𝑠𝑑𝑊_𝑠𝑑𝑡,𝑌_𝑛+1=𝑌_𝑛+𝑔𝑡_𝑛,𝑌_𝑛1−𝑡_𝑛ℎ_𝑛+1+𝜎Δ𝑊_𝑛+1. (19)(20) Label the two integral terms as 𝐿₁ and 𝐿₂, respectively, and let 𝜖_𝑛+1=𝑋_𝑛+1−𝑌_𝑛+1. We have 𝜖_𝑛+1=𝜖_𝑛+ℎ_𝑛+11−𝑡_𝑛⋅𝑔𝑡_𝑛,𝑋_𝑛−𝑔𝑡_𝑛,𝑌_𝑛+𝐿₁+𝐿₂. (21) Using the relation 𝔼𝜖²_𝑛+1=𝔼𝜖²_𝑛+𝔼𝜖_𝑛+1−𝜖_𝑛²+2𝔼𝜖_𝑛+1−𝜖_𝑛𝜖_𝑛 (22) we have 𝔼𝜖²_𝑛+1≤𝔼𝜖²_𝑛𝐶+2⋅²¹ℎ²_𝑛+11−𝑡_𝑛²⋅𝔼𝜖²_𝑛+4𝔼𝐿²₁+4𝔼𝐿²₂𝐶+2⋅¹⋅ℎ_𝑛+11−𝑡_𝑛⋅𝔼𝜖²_𝑛+2𝔼𝜖_𝑛⋅𝐿₁+2𝔼𝜖_𝑛⋅𝐿₂. (23) We have also used the inequality 𝔼(𝑋+𝑌)²≤2𝔼(𝑋²)+2𝔼(𝑌²) and the Lipschitz condition (17) here. Since 𝐿₂ is a martingale with respect to the natural filtration, 𝔼[𝜖_𝑛⋅𝐿₂]=0. Now we look at the second moments of the integral terms 𝐿₁ and 𝐿₂. First we use Cauchy-Schwartz inequality twice on 𝐿₁ and the bounds imposed in (16a)–(16d) 𝔼𝐿²₁=𝔼^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛(𝐴𝑎)_𝑠𝑑𝑠𝑑𝑡²≤ℎ_𝑛+1⋅𝔼^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛(𝐴𝑎)_𝑠𝑑𝑠²𝑑𝑡≤ℎ²_𝑛+1⋅^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛𝔼(𝐴𝑎)²_𝑠𝑑𝑠𝑑𝑡≤ℎ²_𝑛+1⋅^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛𝔼𝑔𝑠,𝑋_𝑠(1−𝑠)²+𝑔₁𝑠,𝑋_𝑠()+𝑔1−𝑠𝑠,𝑋_𝑠𝑔₂𝑠,𝑋_𝑠(1−𝑠)²+𝜎²𝑔₂₂(𝑠,𝑋_𝑠)2(1−𝑠)²𝑑𝑠𝑑𝑡≤ℎ²_𝑛+1⋅^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛𝐷₁(1−𝑠)⁴+𝐷₂(1−𝑠)³+𝐷₃(1−𝑠)²𝑑𝑠𝑑𝑡 (24) for some 𝐷₁, 𝐷₂, and 𝐷₃ that do not depend on 𝑛 or ℎ_𝑛+1. It follows that 𝔼𝐿²₁≤ℎ⁴_𝑛+1⋅𝐷₁1−𝑡_𝑛+1⁴+ℎ⁴_𝑛+1⋅𝐷₂1−𝑡_𝑛+1³+ℎ⁴_𝑛+1⋅𝐷₃1−𝑡_𝑛+1²≤𝐷₄ℎ⁴_𝑛+11−𝑡_𝑛+1⁴, (25) where 𝐷₄ is independent of 𝑛 and ℎ_𝑛+1. Now for 𝔼[𝐿²₂] term we use the Cauchy-Schwartz inequality and the Itô isometry 𝔼𝐿²₂=𝔼^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛(𝐵𝑎)_𝑠𝑑𝑊_𝑠𝑑𝑡²=𝔼^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛𝜎⋅𝑔₂(𝑠,𝑋_𝑠)(1−𝑠)𝑑𝑊_𝑠𝑑𝑡²≤ℎ_𝑛+1⋅^𝑡𝑛+1_𝑡𝑛𝔼^𝑡_𝑡𝑛𝜎⋅𝑔₂(𝑠,𝑋_𝑠)(1−𝑠)𝑑𝑊_𝑠²𝑑𝑡≤ℎ_𝑛+1^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛𝜎𝐶₁(1−𝑠)²𝑑𝑠𝑑𝑡≤𝐷₅ℎ³_𝑛+11−𝑡_𝑛+1², (26) where 𝐷₅ is independent of 𝑛 and ℎ_𝑛+1.

We use Cauchy-Schwartz inequality and the inequality 𝔼𝜖²_𝑛^1/2≤𝑘+𝑘⁻¹𝔼𝜖²_𝑛 (27) with 𝑘=ℎ_𝑛+1/(1−𝑡_𝑛+1) to estimate 𝔼[𝜖_𝑛⋅𝐿₁]: 𝔼𝜖_𝑛⋅𝐿₁≤𝔼𝐿²₁^1/2⋅𝔼𝜖²_𝑛^1/2≤𝐷₄ℎ²_𝑛+11−𝑡_𝑛+1²ℎ_𝑛+11−𝑡_𝑛+1+1−𝑡_𝑛+1ℎ_𝑛+1𝔼𝜖²_𝑛≤𝐷₆ℎ_𝑛+11−𝑡_𝑛+1𝔼𝜖²_𝑛+𝐷₇ℎ³_𝑛+11−𝑡_𝑛+1³. (28)

For simplicity, we assume 𝑌₀=𝑥₀. For fixed step size, ℎ_𝑛=ℎ=1/𝑁, there is no uniform control over the local truncation error since, close to 𝑡=1, 1−𝑡_𝑁−1=ℎ. If we fix 𝛿>0 and only consider 𝑡_𝑛+1<1−𝛿, then we have the following inequality: 𝔼𝜖²_𝑛+1≤1+𝑚₁ℎ𝛿𝔼𝜖²_𝑛+𝑚₂ℎ³𝛿³, (29) where 𝑚₁ and 𝑚₂ are some constants which do not depend on ℎ or 𝛿. Using Lemma 2, we have 𝔼𝜖²_𝑛𝑚≤exp¹𝛿𝑚₂ℎ²𝑚₁𝛿². (30) For fixed 𝛿>0, this is the original case where the coefficients of the SDE have bounded derivatives. However, if 𝛿↓0, this estimate shows no decrease in global error unless ℎ decays exponentially relative to 𝛿.

To produce a better estimate and a more efficient algorithm, we consider variable step sizes defined as follows. First fix ℎ^∗, 0<ℎ^∗<1. Then define step sizes ℎ_𝑛 and node points 𝑡_𝑛 using ℎ^∗:ℎ_𝑛=ℎ^∗1−𝑡_𝑛−1,𝑡_𝑛=𝑡_𝑛−1+ℎ^∗1−𝑡_𝑛−1. (31a)(31b)This implies that𝑡_𝑛=1−1−ℎ^∗𝑛,ℎ_𝑛=ℎ^∗1−ℎ^∗𝑛−1. (32a)(32b)This specific definition allows us to estimate the global error using ℎ^∗, and therefore we can control ℎ^∗ to control the global error. Using these step sizes, we have the error estimates𝔼𝐿²₁≤𝐷₄ℎ^∗4,𝔼𝐿²₂≤𝐷₅ℎ^∗3,𝔼𝜖_𝑛⋅𝐿₁≤𝐷₆ℎ^∗𝔼𝜖²_𝑛+𝐷₇ℎ^∗3, (33a)(33b)(33c)where 𝐷₄, 𝐷₅, 𝐷₆, and 𝐷₇ are positive constants which do not depend on 𝑛 or ℎ^∗. Using these estimates in (23), we have 𝔼𝜖²_𝑛+1≤1+𝑚₁ℎ^∗𝔼𝜖²_𝑛+𝑚₂ℎ^∗3, (34) where 𝑚₁ and 𝑚₂ are constants which do not depend on 𝑛 or ℎ^∗. Using Lemma 2, we have the following.

Theorem 3. Given that the SDE in (15a)-(15b) satisfies the assumptions (16a)–(16d) and 𝔼[(𝑥₀−𝑌₀)²]≤𝐷₀ℎ^∗2, then the above constructed algorithm (20) with variable step sizes (32b) has 𝔼[𝜖²_𝑛]=𝑂(ℎ^∗2), as ℎ^∗↓0, uniformly in 𝑛 for 𝑡_𝑛<1−𝛿<1, and thus the global pointwise error for the above proposed algorithm is of order 𝑂(ℎ^∗).

Proof. If we have 𝑁 steps, (32b) gives (1−ℎ^∗)^𝑁=𝛿, and thus 𝑁=ln𝛿ln1−ℎ^∗≤−ln𝛿ℎ^∗,asℎ^∗↓0. (35) Then by using Lemma 2 on (34) and this 𝑁 we have 𝔼𝜖²_𝑛≤exp−𝑚₁ℎ^∗⋅ln𝛿ℎ^∗⋅𝑚₂ℎ^∗3𝑚₁ℎ^∗+𝐷₀ℎ^∗2−𝑚₂ℎ^∗3𝑚₁ℎ^∗≤𝐷𝛿^𝑚1⋅ℎ^∗2, (36) where 𝐷 and 𝑚₁ are constants that do not depend on 𝑛, ℎ^∗, or 𝛿.

Corollary 4. Given a positive constant 𝑝<1 and the SDE (15a)-(15b) satisfying the assumptions in Theorem 3, the above constructed numerical algorithm has global pointwise error 𝑂(ℎ^∗𝑝), as ℎ^∗↓0, when 𝛿=(ℎ^∗)^{2(1−𝑝)/𝑚1}, where 𝑚₁ is the constant in (34) which only depends on the bounds imposed on 𝑔 and the moments of the process.

Proof. Letting 𝛿=(ℎ^∗)^{2(1−𝑝)/𝑚1}, we have 𝔼𝜖²_𝑛≤𝐷ℎ^∗2𝑝. (37)

With a simple example, we show that we cannot expect a much better result. Consider the SDE: 𝑑𝑋_𝑡=𝑎1−𝑡𝑑𝑡+𝜎𝑑𝑊_𝑡,𝑋₀=0, (38) where 𝑊₀=0. This equation has solution 𝑋_𝑡=−𝑎ln(1−𝑡)+𝜎𝑊_𝑡. (39) The variable step size algorithm generates the approximation 𝑌_𝑛=𝑎𝑛ℎ^∗+𝜎𝑊_𝑡𝑛. (40) It follows then that the maximum value of |𝑋_𝑡𝑛−𝑌_𝑛|=𝑎𝑛(−ln(1−ℎ^∗)−ℎ^∗) occurs at 𝑁 where (1−ℎ^∗)^𝑁=𝛿 and 𝑋_𝑡𝑁−𝑌_𝑁()ℎ=𝑎ln𝛿−1−^∗ln1−ℎ^∗∼𝑎2ℎ^∗ln𝛿. (41) Thus if 𝛿=ℎ^∗𝑞, |𝑋_𝑡𝑁−𝑌_𝑁|∼(𝑞𝑎/2)ℎ^∗lnℎ^∗ as ℎ^∗→0 and is thus 𝑜(ℎ^∗𝑝) for all 𝑝<1.

4. Numerical Approximation of Brownian Bridge with Variable Step Sizes

We apply the above approximation to the Brownian bridge introduced in (3a)-(3b). 𝑔(𝑥)=1−𝑥 satisfies the Lipschitz and linear growth conditions assumed. We also show that we can relax the bounded moments conditions for this specific case. We use the BB starting at 0 and converging to 1 as 𝑡→1. Consider the SDE: 𝑑𝑋_𝑡=1−𝑋_𝑡1−𝑡𝑑𝑡+𝑑𝑊_𝑡,0≤𝑡<1,𝑋₀=0, (42a)(42b)where 𝑊_𝑡 is the standard Brownian motion. Then using (10) we get 𝑋_𝑛+1=𝑋_𝑛+1−𝑋_𝑛1−𝑡_𝑛ℎ_𝑛+1+Δ𝑊_𝑛+1+^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛11−𝑠𝑑𝑊_𝑠𝑑𝑡. (43) Note that (𝐴𝑎)_𝑠 term is zero and that gives us the opportunity to relax the bounded moments condition 𝐴𝑎=𝜕𝑎𝜕𝑡+𝑎𝜕𝑎+1𝜕𝑥2𝑏²𝜕²𝑎𝜕𝑥²=𝑏−𝑥(1−𝑡)²−𝑏−𝑥(1−𝑡)²+12𝑏²⋅0=0. (44) Using (11) we get 𝑌_𝑛+1=𝑌_𝑛+1−𝑌_𝑛1−𝑡_𝑛ℎ_𝑛+1+Δ𝑊_𝑛+1. (45) Now define 𝜖_𝑛+1=𝑋_𝑛+1−𝑌_𝑛+1 to get 𝜖_𝑛+1=𝜖_𝑛−𝜖_𝑛1−𝑡_𝑛ℎ_𝑛+1+^𝑡𝑛+1_𝑡𝑛^𝑡_𝑡𝑛11−𝑠𝑑𝑊_𝑠𝑑𝑡. (46) Proceeding as before we have that the second moment of the integral term is of order ℎ^∗3. Now using (22) and analyzing it in the same way as in the previous section, we get 𝔼𝜖²_𝑛+1≤𝔼𝜖²_𝑛+ℎ^∗2𝔼𝜖²_𝑛−2ℎ^∗𝔼𝜖²_𝑛+𝑂ℎ^∗3. (47) For ℎ^∗<1 we get 𝔼𝜖²_𝑛+1≤1−ℎ^∗𝔼𝜖²_𝑛+𝑂ℎ^∗3 (48) which gives the global error 𝑂(ℎ^∗) using Lemma 2. Note that, in this case, 𝔼𝜖²_𝑛≤𝐷ℎ^∗2 (49) for all 𝑛, independent of 𝛿. However, making 𝛿 small relative to ℎ^∗ increases the number of steps significantly. For example, if we set 𝛿=ℎ^∗2, it takes about twice as many steps to reach 1−𝛿 as in the case 𝛿=ℎ^∗.

For this example, we can verify the error explicitly. Using Algorithm (45), we have 𝑌_𝑛=𝑡_𝑛+^𝑛_𝑘=11−ℎ^{∗𝑛−𝑘}Δ𝑊_𝑘. (50) Using this with the explicit solution (4), we calculate 𝔼𝜖²_𝑛=1−ℎ^∗𝑛1−1−ℎ^∗𝑛×1+2ln1−ℎ^∗ℎ^∗+11−ℎ^∗ (51) and thus the maximum value of 𝔼[𝜖²_𝑛]∼ℎ^∗2/12 as ℎ^∗→0.

Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors would like to thank the referee for insightful comments and help in improving this paper.