Dynamic Systems and Applications 15 (2006) 287-300

APPROXIMATE SOLUTION OF REFLECTING BROWNIAN MOTION USING PENALTY METHOD AND NUMERICAL APPLICATION TO IMPERFECT ELASTIC BARRIER

S. KANAGAWA, Y. SAISHO, AND H. UESU

Department of Mathematics, Faculty of Engineering, Musashi Institute of Technology, Tokyo, 158-8557 Japan Applied Mathematics Research Group, Graduate School of Engineering, Hiroshima University, Hiroshima, 739-8527 Japan Department of Mathematics, Faculty of Education, Waseda University, Tokyo, 169-88050 Japan

ABSTRACT. We investigate the error of the Euler-Maruyama approximate solution of reflecting Brownian motion using the penalty method. The approximate solution is constructed from not only i.i.d. random variables but also dependent sequence of random variables, e.g. Gaussian sequence, mixing sequence, etc. Further we show some numerical applications. AMS (MOS) Subject Classification. 60F17, 60G50, 60H35, 65C30. Key words: Euler-Maruyama approximation, reflecting Brownian motion, stochastic differential equation, Monte Carlo simulation, pseudo-random numbers, imperfect elastic barrier

Dedicated to the memory of Professor Hiroshi Suda

1. INTRODUCTION

In Kanagawa-Saisho [7] we considered a reflecting Brownian motion defined by some stochastic differential equations with boundary conditions on multi-dimensional domains (so-called Skorohod SDE) and obtained the error estimation of the approximate solution of the SDE. The reflecting Brownian motion is applied in many fields. Recently in mathematical finance, not only in science or engineering, but the reflecting Brownian motion is used to estimate the price movement of stocks and bonds. As for the construction of such approximate solution there are two approaches. Tanaka [16], Lions-Sznitman[9], and Saisho[12] constructed Skorohod equations using the projection π on the boundary. The reflecting path x(t) is defined for given function w(t) by the following manner: Define a step function w_{n}(t) by discretization of w(t) and construct the reflecting step function x_{n}(t) for w_{n}(t) as follows. For a domain D ⊂ R^{d}, if x_{n}(t_{0}) ∈ \overline{D} and y_{n}(s) ≡ x_{n}(t_{0}) + w_{n}(s) − w_{n}(t_{0}) ∈ \overline{D}, t_{0} ≤ s ≤ t for some t_{0} and t, put x_{n}(s) = y_{n}(s), t_{0} ≤ s ≤ t. On the other hand, if y_{n}(s) ∈ \overline{D}, t_{0} ≤ s < t and y_{n}(t) ∉ \overline{D}, put x_{n}(t) = π(y_{n}(t)). This type of approximation is often called projection scheme ([12]). A second one is the penalty method type, that is, we approximate Skorohod equations by equations with coefficients of gradient type ([14]). For the purpose of obtaining numerical data by computer simulation of Skorohod SDE’s we need discretized approximate solutions of them. The penalty method is separated into two steps. The first step is to approximate reflecting Brownian motions by usual Brownian motions with a drift term whose coefficients are gradient type in pathwise sense. The second step is to approximate the Brownian motion with such drift term by Euler-Maruyama scheme.

Received April 30, 2004           1056-2176 $15.00 \begin{matrix}c \\ ◯ \end{matrix}Dynamic Publishers, Inc.

In this paper we have two aims. One is to define the approximate solution process using not only a sequence of i.i.d. random variables like increments of Brownian motion but also dependent sequence that does not obey normal distribution. Another one is, to show the advantage of the penalty method, we define a random walk to approximate sample paths of Brownian motion with imperfect elastic barriers, i.e. any path of the Brownian motion does not reflect at the boundary immediately but is absorbed for a short period according to the strength of the path getting out of the boundary. Applying the method we see some examples of numerical computation to simulate such reflecting Brownian motion.

2. SKOROHOD PROBLEM

We define a reflecting Brownian motion satisfying the Skorohod equation,

\begin{cases}x (t) = w (t) + φ (t) , & t ≥ 0 \\ w (0) ∈ D̄. \end{cases}

The Skorohod problem means to find a pair of functions (x, φ) satisfying the Skorohod equation under conditions

• x (t) ∈ \overline{D}, is continuous,

• φ (t) : is continuous and bounded variation, φ (0) = 0,

• φ (t) = ∫_{0}^{t} n (s)d |φ|_{s} , |φ|_{t} = ∫_{0}^{t} I_{{x(s)∈∂D}}(s)d |φ|_{s},

• n (s) ∈ N_{x(s)} if x (s) ∈ ∂D

for given continuous function w. Here, |φ|_{t} is the total variation of φ on [0, t] and

I_{A} (s) = \begin{cases}1, & s ∈ A \\ 0, & s ∉ A. \end{cases}

It is well known that, for the convex domain D, the Skorohod equation can be solved and there exists a reflecting Brownian motion in D. For more details, see [16]. In this paper we call a reflecting Brownian motion when it solves the Skorohod problem. Now, to see that there exist reflecting Brownian motions with non-convex domains, we show some conditions obtained in [12] under which the Skorohod equation can be satisfied. Following [12] we pose two conditions for a domain D ⊂ R^{d}. Define sets of inward normal vectors at x ∈ ∂D by

N_{x,r} = {n : |n| = 1, B (x − rn, r) ∩ D = ∅} , r > 0,

N_{x} = ⋃_{r>0} N_{x,r}, x ∈ ∂D,

where B(y, r) = {z ∈ R^{d} : |z − y| < r}, y ∈ R^{d}, r > 0, and pose two conditions for domains. Condition (A) There exists a constant r_{0} > 0 such that for all x ∈ ∂D

N_{x} = N_{x,r_{0}} ≠ ∅.

Condition (B) There exist constants δ > 0 and β ∈ [1, ∞) such that the following holds: for any x ∈ ∂D there is a unit vector ℓ_{x} with

⟨ℓ_{x}, n⟩ ≥ 1/β for any n ∈ ⋃_{y∈B(x,δ)∩∂D} N_{y},

where ⟨·, ·⟩ is the usual inner product in R^{d}.

In the case when D is convex, Condition (A) is satisfied for any r_{0} > 0. Moreover if D satisfies Condition (A) and dist(x, \overline{D}) < r_{0} there exists a unique π(x) ≡ \overline{x} ∈ \overline{D} with dist(x, \overline{D}) = |x − \overline{x}|. Especially, we have \overline{x} = x for x ∈ D and if x ∉ D we have \overline{x} ∈ ∂D. Existence of the unique solution of the Skorohod equation is guaranteed by the following theorem.

Theorem 1. ([12]) Suppose that D satisfies Conditions (A) and (B). Then for any continuous w, there exists a unique solution (x, φ) of the Skorohod problem.

3. PENALTY METHOD AND EULER-MARUYAMA SCHEME

It is known that if a domain D (⊂ R^{d}) satisfies Condition (A), there is a function U : R^{d} → R_{+} satisfying the following conditions ([14]):

• U ∈ C^{1}(R^{d})

• U(x) = |x − \overline{x}|^{2} if dist(x, \overline{D}) ≤ r_{0}/2

• ∇U is bounded and Lipschitz continuous, i.e. for any x and y,

(3.1)  |∇U (x) − ∇U (y)| ≤ L |x − y| , |∇U (x)| ≤ L.

Now for m > 0 we consider the SDE

(3.2)  X_{m} (t) = x_{0} + B(t) − \frac{m}{2} ∫_{0}^{t} ∇U (X_{m} (s)) ds, 0 ≤ t ≤ T

and the Skorohod equation X(t) = x_{0} + B(t) + Φ(t), X(0) = x_{0} ∈ \overline{D}. Then it is easy to see that the SDE (3.2) has a unique solution. We also have the following results.

Theorem 2. ([14]) Suppose D satisfies Conditions (A) and (B). Then we have as m → ∞

max_{0≤t≤T} |X_{m} (t) − X (t)| → 0 a.s.

Let △_{n} be a partition of the interval [0, T]:

△_{n} : 0 = t_{0} < t_{1} < t_{2} < ⋯ < t_{n} = T

and put

|△_{n}| = max_{1≤k≤n} |t_{k} − t_{k−1}|, △B(t_{k}) = B(t_{k}) − B(t_{k−1}), △t_{k} = t_{k} − t_{k−1}.

We note that {B(t)} is a standard Brownian motion defined on a probability space (Ω, F, P). Then the Euler-Maruyama scheme for the SDE (3.2) is given by for each 1 ≤ k ≤ n

X_{m}^{n} (t) = x_{k−1} + \frac{t − t_{k−1}}{t_{k} − t_{k−1}} (x_{k} − x_{k−1}) , t_{k−1} ≤ t ≤ t_{k},

where

x_{k} = x_{0} + B (t_{k}) − \frac{m}{2} ∑_{i=1}^{k} ∇U (x_{i−1}) Δt_{i}.

This approximation scheme is first given in [10] and convergence rate is investigated in [1], [3] – [6], [8] and [15].

Theorem 3. ([7]) For any ε > p/2, 2 ≤ p < ∞

E [ max_{0≤t≤T} |X_{m}^{n} (t) − X_{m} (t)|^{p}] = o (|Δ_{n}|^{p/2} (− log |Δ_{n}|^{ε})) , n → ∞.

From the above theorem we have the following result.

Theorem 4. ([7]) Let Δ_{n} = O (\frac{1}{n}) as n → ∞. Then we have for any δ, δ′ > 0

n^{\frac{1}{2} −δ} max_{0≤t≤T} |X_{m}^{n} (t) − X_{m} (t)| → 0 a.s., n → ∞,

m^{\frac{1}{4} −δ′} max_{0≤t≤T} |X_{m} (t) − X (t)| → 0 a.s., m → ∞.

4. RANDOM WALKS VIA DEPENDENT RANDOM VARIABLES TO APPROXIMATE X_{m} (t)

Throughout this section and the rest we will assume that a domain D satisfies Conditions (A) and (B). Since X_{m} (t) is defined via a Brownian motion which have to be approximated by random walks when computer simulation is done, we now construct a sequence of random functions {Y_{m}^{n} (t) , n ≥ 1} to approximate X_{m} (t) using not only increments of Brownian motion but also dependent sequence of random variables, e.g. Markov chains, mixing random variables and so on. The following result is a key theorem to estimate the error between Y_{m}^{n} (t) and X_{m} (t). Let {Y_{m}^{n} (t) , 0 ≤ t ≤ T } be a random function defined by semimartingale of sums of a triangular array of random variables {η_{k}^{n}, 1 ≤ k ≤ n}_{n≥1} on a probability space (Ω, F, P) as follows. For each 1 ≤ k ≤ n

Y_{m}^{n} (t) = y_{k−1} + \frac{t − t_{k−1}}{t_{k} − t_{k−1}} (y_{k} − y_{k−1}) , t_{k−1} ≤ t ≤ t_{k},

where

σ_{k}^{2} = V ar \{∑_{i=1}^{k} ξ_{k}^{n}\} , t_{k} = \frac{σ_{k}^{2}}{σ_{n}^{2}} T

and

\begin{cases}y_{0} = x_{0}, \\ y_{k} = x_{0} + ∑_{i=1}^{k} η_{i}^{n} − \frac{m}{2} ∑_{i=1}^{k} ∇U (y_{i−1}) Δt_{k}, & k = 1, ⋯, n. \end{cases}

Theorem 5. For each m > 0, there exist a Brownian motion {B (t), 0 ≤ t ≤ T } and a triangular array of random variables {ξ_{k}^{n}, 1 ≤ k ≤ n}_{n≥1} on the common probability space (Ω, F, P) such that for any p ≥ 2, as n → ∞

E ( max_{0≤t≤T} |Y_{m}^{n} (t) − X_{m} (t)|^{p})

= E ( max_{1≤i≤n} |\frac{σ}{\sqrt{T}} ∑_{j=1}^{i} ξ_{j}^{n} − B (\frac{i}{n})|^{p}) + O (n^{−p/4}) .

We show easily the following result from Theorem 5 and the Borel-Cantelli lemma.

Theorem 6. If there exist a Brownian motion {B (t) , 0 ≤ t ≤ T } and a triangular array of random variables with as n → ∞

(4.1)  E ( max_{1≤i≤n} |\frac{σ}{\sqrt{T}} ∑_{j=1}^{i} ξ_{j}^{n} − B (\frac{i}{n})|^{p}) = O (n^{−p/4}) ,

then, for any m > 0, we have

max_{0≤t≤T} |Y_{m}^{n} (t) − X_{m} (t)| → 0 as n → ∞ a.s.

5. PROOF OF THEOREM 5

For simplicity put ξ_{i}^{n} = ξ_{i}, σ_{k}^{2} = \frac{k}{n} and T = 1. In what follows K’s are different positive constants independent of n in different equations. We first define two random processes {X̄_{n} (t) , 0 ≤ t ≤ 1} and {Ȳ_{n} (t) , 0 ≤ t ≤ 1} as follows. Let {ζ_{k}, 1 ≤ k ≤ M} and {η_{k}, 1 ≤ k ≤ M} be random variables defined by

ζ_{k} := ∑_{i=(k−1)[n^{1/2}]+1}^{k[n^{1/2}]} \frac{ξ_{i}}{\sqrt{n}}, 1 ≤ k ≤ M − 1,

ζ_{M} := ∑_{i=(M−1)[n^{1/2}]+1}^{n} \frac{ξ_{i}}{\sqrt{n}},

η_{k} := B (k [n^{1/2}]/n) − B (((k − 1) [n^{1/2}] + 1)/n) , 1 ≤ k ≤ M − 1,

and

η_{M} := B (1) − B (((M − 1) [n^{1/2}] + 1)/n) ,

where M := [n/[n^{1/2}]]+1. Define {X̄_{n} (t) , 0 ≤ t ≤ 1} and {Ȳ_{n} (t) , 0 ≤ t ≤ 1} by

\begin{cases}X̄_{n} (t) := u_{k}, ((k − 1) [n^{1/2}] + 1)/n ≤ t < k [n^{1/2}]/n, 1 ≤ k ≤ M − 1 \\ X̄_{n} (1) := u_{M}, \end{cases}

\begin{cases}Ȳ_{n} (t) := v_{k}, ((k − 1) [n^{1/2}] + 1)/n ≤ t < k [n^{1/2}]/n, 1 ≤ k ≤ M − 1 \\ Ȳ_{n} (1) := v_{M}, \end{cases}

where for 1 ≤ k ≤ M

u_{k} := X_{0} + η_{j} + ∑_{j=1}^{k} ∇U (u_{j−1}) [n^{1/2}]/n,

v_{k} := X_{0} + ζ_{j} + ∑_{j=1}^{k} ∇U (v_{j−1}) [n^{1/2}]/n.

From the definitions of X_{n} (t) and Y_{n} (t) (5.1)

E ( max_{0≤u≤t} |X_{m} (u) − Y_{m}^{n} (u)|^{p})^{1/p} ≤ E ( max_{0≤u≤t} |X_{m} (u) − X̄_{m}^{n} (u)|^{p})^{1/p}

+E ( max_{0≤u≤t} |X̄_{m}^{n} (u) − Ȳ_{m}^{n} (u)|^{p})^{1/p} + E ( max_{0≤u≤t} |Ȳ_{m}^{n} (u) − Y_{m}^{n} (u) |^{p})^{1/p}

:= I_{1}^{1/p} + I_{2}^{1/p} + I_{3}^{1/p}.

We first estimate I_{2}. For \frac{k[n^{1/2}]}{n} ≤ t ≤ \frac{(k+1)[n^{1/2}]}{n}, 1 ≤ k ≤ M, we have from (1.3)

(5.2) I_{2} ≤ E( max_{1≤i≤k} |u_{i} − v_{i}|^{p})

≤ KE ( max_{1≤i≤k} |∑_{j=1}^{i} η_{j} − ∑_{j=1}^{i} ζ_{j}|^{p})

+ KE ( max_{1≤i≤k} |∑_{j=1}^{i} (∇U (u_{j−1}) − ∇U (v_{j−1})) \frac{[n^{1/2}]}{n}|^{p})

=: I_{21} + I_{22}.

On the other hand by Doob’s inequality and (3.1) it is easy to see that

(5.3)  I_{22} ≤ K ∫_{0}^{t} E ( max_{0≤u≤s} |X̄_{m}^{n} (u) − Ȳ_{m}^{n} (u) |^{p}) ds.

From (4.2) and (4.3), for any 0 ≤ t ≤ 1

E ( max_{0≤u≤t} |X̄_{m}^{n} (u) − Ȳ_{m}^{n} (u)|^{p}) ≤ E ( max_{1≤i≤k} |∑_{j=1}^{i} η_{j} − ∑_{j=1}^{i} ζ_{j}|^{p})

+K ∫_{0}^{t} E ( max_{0≤u≤s} |X̄_{m}^{n} (u) − Ȳ_{m}^{n} (u)|^{p}) ds.

Thus by Gronwall’s lemma

(5.4)  I_{2} ≤ KE ( max_{1≤i≤k} |∑_{j=1}^{i} η_{j} − ∑_{j=1}^{i} ζ_{j}|^{p}) .

As for I_{1} and I_{3}, from Lemmas 5 and 6 in Ikeda-Watanabe (1981), we have as n → ∞

(5.5)  I_{1} = O (n^{−p/4}) and I_{3} = O (n^{−p/4}) .

From (5.1), (5.4) and (5.5) we conclude the proof of the theorem.

6. EXAMPLES OF Y_{m}^{n} (t)

There are many kinds of sequences of weakly dependent random variables satisfying (4.1), e.g. mixing random variables, Gaussian sequence of random variables and so on. Of course sums of these sequences are not semimartingales, however they can be approximated via martingale differences using the Skorohod embedding theorem. The detail of such approximation is shown in Philipp-Stout [11]. We interest in the case of dependent random variables which define Y_{m}^{n} (t) because of not only the mathematical aspect but also some applications for it. For example, in the case of financial markets, many trades of bonds or stocks are not done independently each other. Usually a trader buys or sells them according to other traders’ behavior. Now, according to [11], we show some examples of dependent random variables for which Theorem 6 holds.

Example 1. Stationary ϕ-Mixing Sequence Let {ξ_{k}^{n}, 1 ≤ k ≤ n}_{n≥1} be a triangular array of strictly stationary ϕ-mixing random variables with the ϕ-mixing coefficients {ϕ_{n} (k)}_{n≥1} satisfying

|P (A ∩ B) − P (A) P (B)| ≤ ϕ_{n} (k) P (A) ,

for A ∈ F {ξ_{1}^{n}, ⋯, ξ_{m}^{n}} and B ∈ F {ξ_{m+k}^{n}, ⋯, ξ_{n}^{n}} for each n uniformly in m. Suppose that

lim_{n→∞} ∑_{k=1}^{n} ϕ_{n} (k) < ∞,

Then we have

σ^{2} := lim_{n→∞} E [(∑_{k=1}^{n} ξ_{k}^{n})^{2}] < ∞.

When σ^{2} > 0, {ξ_{k}^{n}} satisfies (4.1).

As for limit theorems for other types of mixing sequences, see e.g. [17].

Example 2. Gaussian Sequence Let {ξ_{k}^{n}, 1 ≤ k ≤ n}_{n≥1} be a triangular array of Gaussian random variables satisfying

E [( ∑_{k=m+1}^{m+n} ξ_{k}^{n})^{2}] = σ^{2} + O (n^{−ε}) , n → ∞,

uniformly in m for some ε > 0. Let σ^{2} = 1, then {ξ_{k}^{n}} satisfies (4.1).

Example 3. Lacunary Trigonometric Series Let {n_{k}, k ≥ 1} be a lacunary sequence of positive integers satisfying

\frac{n_{k+1}}{n_{k}} ≥ q > 1, k ≥ 1,

for some q > 1. Suppose that the sequence

{cos 2πn_{k}ω, k ≥ 1}

is defined on a probability space {[0, 1) , B, P}, where P is a Lebesgue measure. Put

ξ_{k}^{n} = \frac{1}{σ_{n}} cos 2πn_{k}ω, 1 ≤ k ≤ n, σ_{n}^{2} = V ar \{∑_{k=1}^{n} cos 2πn_{k}ω\} .

Then {ξ_{k}^{n}} satisfies (4.1).

7. NUMERICAL APPLICATIONS

In this section, to show the advantage of the penalty method, we consider sample paths of some reflecting Brownian motion with imperfect elastic barriers, i.e. any path of the Brownian motion is absorbed for a short period according to the strength of the path getting out of the boundary instead of reflecting at the boundary immediately. Such situation is observed when some electromagnetic wave reaches the ionosphere of the earth. This problem was suggested by Professor I. Nagano at Kanazawa University who researches, with NASA, the analysis of wave data observed by scientific satellites, Akebono (EXOS-D) and GEOTAIL in near-Earth space plasmas. As for the hard boundary case we show in [7] a reflecting Brownian motion with a hard non-convex domain whose shape is a concentric circle like the atmosphere of the earth. Furthermore as for an application to interacting particle system of hard balls, see [13].

Now we first treat a reflecting Brownian motion defined on a concentric circle type domain with imperfect elastic barriers.

Example 4. D = {x = (x, y) : 16 < x^{2} + y^{2} < 25} ⊂ R^{2}. The Brownian motion starts from (4, 0) ∈ R^{2} with different m’s which control the hardness of the boundary. The following figures represent paths of X_{m}^{n} starting at (4, 0) ∈ R^{2} with m = 10, 000, n = 2, 500 until 10,000 steps (Fig. 5) and m = 2, 500, n = 2, 500 (Fig. 6), m = 250, n = 2, 500 (Fig. 7) and m = 5, n = 2, 500 (Fig. 8), respectively.

Fig. 5  Fig. 6  Fig. 7  Fig. 8

Since each m in the above example controls the power to pull the Brownian motion path into the domain D from the outside of it, the length of the period in which the path is absorbed in the outside of D depends on m. We next observe the relation between m and the proportion of the numbers of X_{m}^{n} (t_{k}) located outside of D to n steps.

Throughout the following examples {η_{k}^{n}, 1 ≤ k ≤ n} is a sequence of normally distributed random variables with η_{k}^{n} ∼ N (0, (\frac{1}{50})^{2}). The domain of the first case is the same as Example 4. However the initial point and the distribution of η_{k}^{n} are different in these cases, respectively.

Example 5. D = {x = (x, y) : 16 < x^{2} + y^{2} < 25} ⊂ R^{2}. The Brownian motion starts from (0, 4.5) ∈ R^{2} with m = 0.01, 0.1, 0.5, 1, 2. η_{k}^{n} ∼ N (0, (\frac{1}{50})^{2}) implies that Δ_{k} = \frac{1}{50} for each 1 ≤ k ≤ n. Therefore, putting

n = 10000, {X_{m}^{10000} (t) , 0 ≤ t ≤ 200} approximates {X_{m} (t) , 0 ≤ t ≤ 200}. The following figure represents a path of X_{m}^{n} starting at (0, 4.5) ∈ R^{2} with m = 0.1, n = 10, 000 (Fig. 9). At Fig. 9 we observe that the proportion of X_{0.1}^{10,000} (t_{k}) located outside of D to 10, 000 steps (n = 10, 000) is 0.0511 which means the number of steps outside of the domain D is 511.

Fig. 9

The next table is a relation between m = 0.01, 0.1, 0.5, 1, 2 and the averages of proportions of X_{m}^{10,000} (t_{k}) outside of the domain D to 10, 000 steps for 100 sample paths, respectively.

m average of proportions
0.01  0.20
0.1  0.022
0.5  0.029
1  0.017
2  0.0090

Example 6. The domain D = {x = (x, y) : x^{2} + y^{2} < 25} is the inside of a circle. The following figure is the Brownian motion which starts from (0, 5) ∈ R^{2} with m = 0.5, n = 10, 000 (Fig. 10).

Fig. 10

The next table is a relation between m = 0.01, 0.1, 1, 2 and the averages of proportions of X_{m}^{10,000} (t_{k}) located outside of the domain D to 10, 000 steps for 100 sample paths, respectively.

m average of proportions
0.01  0.096
0.1  0.032
1  0.0081
2  0.0037

Example 7. The domain is D = \{(x, y) ; 0 < y < 3 − \sqrt{10 − (x − 3)^{2}},

0 < x < 2}. The following figure is the Brownian motion which starts from (0.2, 1.5) ∈ D with m = 0.5, n = 1, 000 (Fig. 11).

Fig.11

The next table is a relation between m = 0.1, 0.5, 1 and the averages of proportions of X_{m}^{1,000} (t_{k}) located outside of the domain D to 1, 000 steps for 100 sample paths, respectively.

m average of proportions
0.1  0.24
0.5  0.074
1  0.039

ACKNOWLEDGEMENT

The authors would like to express their thanks to Professor I. Nagano for helpful comments about the ionosphere of the earth from which they found an application to some models of a reflecting Brownian motion with imperfect elastic barriers. They are thankful to Professors S. Sathananthan and J. Golec to give the opportunity to present this paper. They would like to express their thanks to the referee and for helpful comments. They are also thankful to Professor K. Inoue for his helpful comments about computer programming.

Research supported in part by Grant-in-Aid Scientific Research (No. 16540124), Ministry of Education, Science and Culture.

REFERENCES

[1] O. Faure, Numerical pathwise approximation of stochastic differential equations, Prépublication of LAMM, Ecole Nationale des Ponts et Chaussees, Paris (1990).

[2] N. Ikeda and S. Watanabe, Stochastic Differential Equations and Diffusion Processes, North-Holland Publ. Co., Amsterdam, Oxford, New York; Kodansha Ltd., Tokyo [3] S. Kanagawa, On the rate of convergence for Maruyama’s approximate solutions of

stochastic differential equations, Yokohama Math. J., 36 (1988), 79–85.

[4] S. Kanagawa, The rate of convergence for approximate solutions of stochastic differential equations, Tokyo J. Math., 12 (1989), 33–48.

[5] S. Kanagawa, Confidence intervals of discretized Euler-Maruyama approximate solutions of SDE’s, Theory, Methods and Applications (Proc. 2nd World Congress of Nonlinear Analysts, Greece, 1996), 30 No.7 (1997), 4101–4104.

[6] S. Kanagawa, Some Remarks on Strong Approximation of Reflecting Brownian Motion Using Penalty Method, Proceedings of Neural, Parallel & Scientific Computations, 2 (2002), 63–70, Atlanta, U.S.A.

[7] S. Kanagawa and Y. Saisho, Strong Approximation of Reflecting Brownian Motion Using Penalty Method and its Application to Computer Simulation, Monte Carlo Methods Appl., 6 (2000), 105–114.

[8] S. Kanagawa and S. Ogawa, Numerical solutions of stochastic differential equations and their applications, to appear in Sugaku Expositions, A.M.S., 18 (2005).

[9] P.L. Lions and A.S. Sznitman, Stochastic differential equations with reflecting boundary conditions, Comm. Pure Appl. Math., 37 (1981), 511–537.

[10] G. Maruyama, Continuous Markov processes and stochastic equations, Rend Circ. Mat. Palermo, 4 (1955), 48–90.

[11] W. Philipp and W. Stout, Almost Sure Invariance Principles for Weakly Dependent Random Variables, Memoir A.M.S., 161 (1975).

[12] Y. Saisho, Stochastic differential equations for multi-dimensional domain with reflecting boundary, Probab. Theory Rel. Fields, 74 (1987), 455–477.

[13] Y. Saisho and H. Tanaka, Stochastic differential equations for mutually reflecting Brownian balls, Osaka J. Math. 23 (1986), 725–740.

[14] Y. Saisho and H. Tanaka, On the symmetry of a reflecting Brownian motion defined by Skorohod’s equation for a multi-dimensional domain, Tokyo J.Math. 10 (1987), 419–435.

[15] H. Schurz, Numerical analysis of SDEs without tears, Handbook of Stochastic Analysis and Applications, D. Kannan and V. Lakshmikantham, eds., 237–359, Marcel Dekker, Basel, 2002.

[16] H. Tanaka, Stochastic differential equations with reflecting boundary condition in convex regions, Hiroshima Math. J., 9 (1979), 163–177.

[17] K. Yoshihara, Weakly Dependent Stochastic Sequences and Their Applications, Vol. 4, Curve Estimation based on Weakly Dependent Data, Sanseido Co. Ltd., Tokyo,