Sim.

DiProc: A Package for Simulation of Diusion

Processes in R
Kamal Boukhetala
Department of Probabilities & Statistics
University of Science and Technology Houari Boumediene
E-mail: kboukhetala@usthb.dz
Arsalane Guidoum
Department of Probabilities & Statistics
University of Science and Technology Houari Boumediene
E-mail: starsalane@gmail.com
Preprint submitted to Journal of Statistical Software (JSS), http://www.jstatsoft.
org/, 25 May 2011.
Abstract
The Sim.DiProc package provides a simulation of diusion processes and the
dierences methods of simulation of solutions for stochastic dierential equations
(SDEs) of the Itos type, in nancial and actuarial modeling and other areas of appli-
cations, for example the stochastic modeling and simulation of pollutant dispersion
in shallow water using the attractive center, and the model of two diusions in at-
traction, which can modeling the behavior of two insects, one attracts the other. The
simulation of the processes of diusion, through stochastic dierential equations to
allow simulated a random variable
c
rst passage time of the particle through a
sphere of radius c, two methods are used in the estimation problem of the probabil-
ity density function of a random variable
c
the histograms and the kernel methods.
The R package Sim.DiProc is introduced, providing a simulation and estimation for
the stationary distribution of the stochastic process that describes the equilibrium of
some dynamics.
Key words : attractive model, diusion process, simulations, stochastic dierential
equation, stochastic modeling, R language.
1
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

1 Introduction
Stochastic dierential equations (SDEs) are a natural choice to model the time evolution
of dynamic systems. These equations have a variety of applications in many disciplines
and can be a powerful tool for the modeling and the description of many phenomena.
Examples of these applications are physics, astronomy, economics, nancial mathematics,
geology, genetic analysis, ecology, neurology, biology, biomedical sciences, epidemiology,
political analysis and social processes, and many other elds of science and engineering.
The stochastic dierential equations, with slight notational variations, are standard in
many books with applications in dierent elds see [23, 15, 1, 25, 14] to name only a few.
The rst is to recall the theory and implement methods for the simulation of paths of
diusion processes {X
t
, t 0} solutions to stochastic dierential equations (SDEs), the
sense that only SDEs with Gaussian noise B
t
are considered i.e., processes for which the
writing
dX
t
dt
= (, t, X
t
) + (, t, X
t
)B
t
With the gaussian noise B
t
is the formal derivative of the standard Wiener process W
t
(i.e.,
B
t
=
dW
t
dt
), the write formally for the process {X
t
, t 0} is
dX
t
= (, t, X
t
)dt + (, t, X
t
)dW
t
,
with some initial condition X
0
.
We seek to motivate further interest in this specic eld by introducing the Sim.DiProc
package [9] (Simulation of Diusion Processes) for use with the statistical programming
environment R [22]; freely available on the Comprehensive R Archive Network (CRAN) at
http://CRAN.R-project.org/package=Sim.DiffProc. There exist a graphical user inter-
face (GUI) for some functions for Sim.DiProc package, the R package Sim.DiProcGUI [10]
a freely available at http://CRAN.R-project.org/package=Sim.DiffProcGUI.
This work is organized as follows. Section 2 gives a simulation for some very popu-
lar stochastic dierential equations models (trajectory of the Brownian motion, Ornstein-
Uhlenbeck process), and present the dierent numerical methods (schemes) of simulation
of SDEs (See [21]), simulation the random variable X
t
at time t by a simulated diusion
processes used the numerical methods and estimate the stationary distribution for X
t
by
histograms and kernel methods. Section 3 the diusion processes are used to modeling
the behavior of the dispersal phenomenon and dynamic models for insects orientation, two
important models can be indicated. The graphical user interface (GUI) for Sim.DiProc
package are provided in Section 4.
2 Diusion processes
We consider the model as the parametric It o stochastic dierential equation
dX
t
= (, t, X
t
)dt + (, t, X
t
)dW
t
, t 0 , X
0
= (1)
2
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

Where {W
t
, t 0} is a standard Wiener process, : [0, T] R R, called the drift
coecient, and : [0, T] R R
+
, called the diusion coecient, are known functions
except the unknown parameters and , R, R and E(
2
) < , the estimation
problems for the parameters and can be seen in some papers [2, 4, 11]. The drift
coecient is also called the trend coecient or damping coecient or translation coecient.
The diusion coecient is also called volatility coecient. Under global Lipschitz and the
linear growth conditions on the coecients and , there exists a unique strong solution of
the above It o SDEs, called the diusion process or simply a diusion, which is a continuous
strong Markov semimartingale. The drift and the diusion coecients are respectively the
instantaneous mean and instantaneous standard deviation of the process. All over the text,
the stochastic dierential equation (1), they are supposed to be measurable.
Assumption 1: (Global Lipschitz) For all x, y R and t [0, T], there exists a
constant K < such that
|(t, x) (t, y)| +|(t, x) (t, y)| < K|x y|.
Assumption 2: (Linear growth) For all x, y R and t [0, T], there exists a
constant C < such that
|(t, x)| +|(t, x)| < C(1 +|x|),
The linear growth condition controls the behaviour of the solution so that X
t
does not
explode in a nite time.
2.1 Simulation of Stochastic Dierential Equations SDEs
In this section, we present some of the well-known and widely used diusion process solu-
tions to the stochastic dierential equation. There are two main objectives in the simulation
of the trajectory of a process solution of a stochastic dierential equation, either interest is
in the whole trajectory or in the expected value of some functional of the diusion process
(moments, distributions, etc) which usually are not available in explicit analytical form.
Numerical Methods are usually based on discrete approximations of the continuous solu-
tion to a stochastic dierential equation. The following examples use the Sim.DiProc
package for simulation some diusion process.
To install the Sim.DiProc package on your version of R, type the following line in the R
console.
R> install.packages("Sim.DiffProc")
If you dont have enough privileges to install software on your machine or account, you
will need the help of your system administrator. Once the package has been installed, you
can actually use it by loading the code with
R> library(Sim.DiffProc)
A short list of help topics, corresponding to most of the commands in the package, is
available by typing
R> library(help = "Sim.DiffProc")
3
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

2.1.1 Simulation of the trajectory of the Brownian motion
Simulation of the trajectory of the Brownian motion: The very basic ingredient of a model
describing stochastic evolution is the so-called Brownian motion or Wiener process. There
are several alternative ways to characterize and dene the Wiener process W = {W
t
, t
0}, and one is the following: it is a Gaussian process with continuous paths and with
independent increments such that W(0) = 0 with probability 1, E(W
t
) = 0, and var(W
t

W
s
) = t s for all 0 s < t. In practice, what is relevant for our purposes is that
W(t) W(s) N(0, t s). Given a xed time increment t > 0, one can easily simulate
a trajectory of the Wiener process in the time interval [t
0
, T]. Indeed, for W
t
it holds true
that
W(t) = W(t) W(0) N(0, t)

t N(0, 1)
and the same is also true for any other increment W(t + t) W(t); i.e.,
W(t + t) W(t) N(0, t)

t N(0, 1)
For i = 0, 1, . . . , N 1, with initial deterministic value x
0
. Usually the time increment
t = t
i+1
t
i
is taken to be constant (i.e.,t = (T t
0
)/N)
R> BMN(N = 10000, t0 = 0, T = 1, C = 1)
R> BMN2D(N = 10000, t0 = 0, T = 1, x0 = 0, y0 = 0, Sigma = 1)
Figure 1: Simulation examples to illustrate the trajectory of the Brownian motion used
the function BMN, and 2-dimensional Brownian motion used BMN2D.
4
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

2.1.2 Simulation a Ornstein-Uhlenbeck or Vasicek process
The Ornstein-Uhlenbeck or Vasicek process is the unique solution to the following stochastic
dierential equation
dX
t
= r( X
t
)dt + dW
t
, X
0
= x
0
, (2)
where is interpreted as the volatility, is the long-run equilibrium value of the process,
and r is the speed of reversion. As an application of the Ito lemma, we can show the explicit
solution of (2) by choosing f(t, x) = xe
rt
, we obtain
X
t
= + (x
0
)e
rt
+
_
t
t
0
e
r(ts)
dW
s
It can be seen that for = 0 the trajectory of X
t
is essentially a negative exponential
perturbed by the stochastic integral. One way of simulating trajectories of the Ornstein-
Uhlenbeck process is indeed via the simulation of the stochastic integral.
R> HWV(N = 1000, t0 = 0, T = 10, x0 = -5, theta = 0, r = 1, sigma = 0.5)
R> HWVF(N = 1000, M = 100, t0 = 0, T = 10, x0 = -5, theta = 0, r = 1,
sigma = 0.5)
Figure 2: Simulated path of the Ornstein-Uhlenbeck process dX
t
= X
t
dt + 0.5dW
t
used
the function HWV, and 100 trajectories of the dX
t
= X
t
dt + 0.5dW
t
used HWVF.
We sketch some very popular SDEs models and example of simulations :
5
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

Models Expressions R code
Arithmetic Brownian Motion dX
t
= dt + dW
t
ABM
Geometric Brownian Motion dX
t
= X
t
dt + X
t
dW
t
GBM
Cox-Ingersoll-Ross dX
t
= (r X
t
)dt +

X
t
dW
t
CIR
Constant Elasticity of Variance dX
t
= X
t
dt + X

t
dW
t
CEV
Radial Ornstein-Uhlenbeck dX
t
= (X
1
t
X
t
)dt + dW
t
ROU
Chan-Karloyi-Logsta-Sanders dX
t
= (r + X
t
)dt + X

t
dW
t
CKLS
Hyperbolic Diusion dX
t
= X
t
(1 + X
2
t
)
0.5
dt + dW
t
Hyproc
Jacobi diusion dX
t
= (X
t
0.5)dt +
_
X
t
(1 X
t
)dW
t
JDP
Table 1: Some parametric families of SDEs models.
2.1.3 Example of use
In particular, it is possible to generate M independent trajectories of the same process
with one single call of the function by just specifying a value for M
R> ABM(N = 1000, t0 = 0, T = 1, x0 = 0, theta = 3, sigma = 2)
R> GBM(N = 1000, T = 1, t0 = 0, x0 = 1, theta = 2, sigma = 0.5)
R> CIR(N = 1000, M = 1, t0 = 0, T = 1, x0 = 1, theta = 0.2, r = 1,
sigma = 0.5)
R> CEV(N = 1000, M = 1, t0 = 0, T = 1, x0 = 1, mu = 0.3, sigma = 2,
gamma = 1.2)
R> ROU(N = 1000, M = 1, T = 1, t0 = 0, x0 = 1, theta = 0.05)
R> CKLS(N = 1000, M = 1, T = 1, t0 = 0, x0 = 1, r = 0.3, theta = 0.01,
sigma = 0.1, gamma = 0.2)
R> Hyproc(N = 1000, M = 1, T = 100, t0 = 0, x0 = 3, theta = 2)
R> JDP(N = 1000, M = 1, T = 100, t0 = 0, x0 = 0, theta = 0.05)
2.2 Numerical Methods for SDEs
The idea is the following given an Ito process {X
t
, 0 t < T} solution of the stochastic
dierential equation
dX
t
= f(X
t
)dt + g(X
t
)dW
t
, X
0
= x
0
, (3)
where W
t
represents the standard Wiener process and initial value x
0
is a xed value. In
many literatures, whose partial list can be seen in the references of the present paper,
numerical schemes for SDE (3) were proposed, which recursively compute sample paths
(trajectories) of solution X
t
at step-points. Numerical experiments for these schemes can
be seen in some papers [16, 20, 24].
In the following, we present numerical schemes. They adopt an equidistant discretiza-
tion of the time interval [t
0
, T] with stepsize
t =
(T t
0
)
N
, for xed natural number N.
6
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

Furthermore,
t
n
= nt, n {1, 2, . . . , N}.
denotes the n-th step-point. We abbreviate X
n
= X
t
n
.
The following three random variables will be used in the (n + 1) time step of the
schemes:
W
n
= W
t
n+1
W
t
n
,
Z
n
=
_
t
n+1
t
n
_
s
t
n
dW
r
ds,
Z
n
=
_
t
n+1
t
n
_
s
t
n
drdW
s
.
They are obtained as sample values of normal random variables using the transforma-
tion [24]
W
n
=
n,1
(t)
1/2
,
Z
n
=
1
2
_

n,1
+

n,2

3
_
(t)
3/2
,
Z
n
=
1
2
_

n,1

n,2

3
_
(t)
3/2
.
and, together with them, we further use

W
n
=
n,2
(t)
1/2
, where
n,1
,
n,2
are mutually
independent N(0, 1) random variables.
2.2.1 Numerical schemes
Euler-Maruyama scheme (Maruyama 1955):
X
n+1
= X
n
+ f
n
t + g
n
W
n
(4)
Milstein scheme (Milstein 1974):
X
n+1
= X
n
+ f
n
t + g
n
W
n
+
1
2
g

n
g
n
((W
n
)
2
t) (5)
Milstein Second scheme (Milstein 1974):
X
n+1
= X
n
+ f
n
t + g
n
W
n
+
1
2
g

n
g
n
((W
n
)
2
t) + f

n
g
n
Z
n
+
_
g

n
f
n
+
1
2
g

n
g
2
n
_
Z
n
+
1
6
(g
2
n
g
n
+ g

n
g
2
n
)((W
n
)
3
3tW
n
)
(6)
7
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

Taylor scheme [20]:
X
n+1
= X
n
+ f
n
t + g
n
W
n
+
1
2
g

n
g
n
((W
n
)
2
t) + f

n
g
n
Z
n
+
_
g

n
f
n
+
1
2
g

n
g
2
n
_
Z
n
+
1
6
(g
2
n
g
n
+ g

n
g
2
n
)((W
n
)
3
3tW
n
)
+
1
2
_
f

n
f
n
+
1
2
f

n
g
2
n
_
(t)
2
(7)
Heun scheme (McShane 1974):
X
n+1
= X
n
+
1
2
[F
1
+ F
2
]t +
1
2
[G
1
+ G
2
]W
n
, (8)
where
F
1
= F(X
n
), G
1
= g(X
n
),
F
2
= F(X
n
+ F
1
t + G
1
W
n
), G
2
= g(X
n
+ F
1
t + G
1
W
n
),
F
x
=
_
f
1
2
g

g
_
(x).
Improved 3-stage Runge-Kutta scheme [24]:
X
n+1
= X
n
+
1
4
[F
1
+ 3F
3
]t +
1
4
[G
1
+ 3G
3
]W
n
+
1
2

3
_
f

n
g
n
g

n
f
n

1
2
g

n
g
2
n
_
t

W
n
,
(9)
where
F
1
= F(X
n
), G
1
= g(X
n
),
F
2
= F
_
X
n
+
1
3
F
1
t +
1
3
G
1
W
n
_
, G
2
= g
_
X
n
+
1
3
F
1
t +
1
3
G
1
W
n
_
,
F
3
= F
_
X
n
+
2
3
F
2
t +
2
3
G
2
W
n
_
, G
3
= g
_
X
n
+
2
3
F
2
t +
2
3
G
2
W
n
_
,
F
x
=
_
f
1
2
g

g
_
(x).
2.2.2 Example of use
The following examples for dierent methods of simulation of SDEs use the snssde function
snssde(N, M, T = 1, t0, x0, Dt, drift, diffusion, Output = FALSE,
Methods = c("SchEuler", "SchMilstein", "SchMilsteinS",
"SchTaylor", "SchHeun", "SchRK3"))
8
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

with : SchEuler (4), SchMilstein (5), SchMilsteinS (6), SchTaylor (7), SchHeun (8),
SchRK3 (9).
Consider for example the stochastic process {X
t
, t 0} solution of
dX
t
= (tX
t
X
3
t
)dt + dW
t
, X
0
= x
0
,
For this process, f(X
t
) = (tX
t
X
3
t
) and g(X
t
) = . Suppose we x an initial value
X
0
= 0 and the set of parameters = 0.03 and = 0.1. The following algorithm can be
used to simulate one trajectory of the process X
t
using the Euler algorithm (gure 3), and
100 trajectories used Milstein scheme (gure 4)
R> f_x <- expression( 0.03 * t * x - x^3 )
R> g_x <- expression( 0.1 )
R> snssde(N = 1000, M = 1, t0 = 0, x0 = 0, Dt = 0.1, drift = f_x,
diffusion = g_x, Methods = "SchEuler")
R> output <- data.frame(time,X)
R> output
Figure 3: Simulated one trajectory of the
process dX
t
= (0.03tX
t
X
3
t
)dt + 0.1dW
t
used Euler algorithm.
Figure 4: Simulated 100 trajectories of the
process dX
t
= (0.03tX
t
X
3
t
)dt + 0.1dW
t
used Milstein scheme.
1 0.0 0.000000e+00
2 0.1 3.162278e-02
3 0.2 6.324555e-06
9
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

. . .
. . .
. . .
1000 99.9 -1.716814e+00
1001 100.0 -1.756944e+00
R> snssde(N = 1000, M = 100, t0 = 0, x0 = 0, Dt = 0.1, drift = f_x,
diffusion = g_x, Methods = "SchMilstein")
2.3 Stationary distribution in the SDEs models
The stationary distribution of the stochastic process that describes the equilibrium of some
dynamics. For these disciplines, the interest is in the shape of the stationary distributions
and the statistical indexes related to them (mean, mode, etc.). The stochastic dierential
equations usually have a linear drift of the form f(x) = r( x) with r > 0. For this
reason, the models are called linear feedback models. The diusion coecient g(x) may be
constant, linearly depending on x, or of polynomial type, leading, respectively, to Gaussian,
Gamma, or Beta stationary distributions.
1
2.3.1 Type N model (Normal Distribution)
For this model, we have:
Drift coecient: f(X
t
) = r( X
t
), r > 0.
Diusion coecient: g(X
t
) =

2, > 0.
Sde: dX
t
= r( X
t
)dt +

dW
t
.
Stationary density: (x) =
1

2
exp
_

(x )
2
2
_
, =

r
.
Statistics: mean, mode = , variance = .
2.3.2 Type G model (Gamma Distribution)
For this model, we have:
1
These correspond to the following Pearson family of distributions: Gamma = type III, Beta = type
I, and Gaussian = limit of type I, III, IV , V , or V I.
10
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

Drift coecient: f(X
t
) = r( X
t
), r > 0.
Diusion coecient: g(X
t
) =
_
X
t
, > 0.
Sde: dX
t
= r( X
t
)dt +
_
X
t
dW
t
.
Stationary density: (x) =
_
x

_
1+

_, =

2r
.
Statistics: mean = , mode = , variance = .
2.3.3 Type B model (Beta Distribution)
For this model, we have:
Drift coecient: f(X
t
) = r( X
t
), r > 0.
Diusion coecient: g(X
t
) =
_
X
t
(1 X
t
), > 0.
Sde: dX
t
= r( X
t
)dt +
_
X
t
(1 X
t
)dW
t
.
Stationary density: (x) =

_
1

_
1

_x
1+

(1 x)
1+
1

, =

2r
.
Statistics: mean = , mode =

1 2
, variance =
(1 )
1 +
.
2.3.4 Example of use
simulation M-sample (M = 100) for the random variable X
t
at time t by a simulated
diusion processes, using the function AnaSimX.
AnaSimX(N, M, t0, Dt, T = 1, X0, v, drift, diff, Output = FALSE,
Methods = c("Euler", "Milstein", "MilsteinS", "Ito-Taylor",
"Heun", "RK3"))
Consider for example the stochastic process the type N : dX
t
= 2(1 X
t
)dt + dW
t
.
R> r = 2
R> theta = 1
R> sigma = 1
R> f_x <- expression(r * (theta - x))
R> g_x <- expression(sqrt(sigma))
R> AnaSimX(N = 1000, M = 100, t0 = 0, Dt = 0.01, T = 10, X0 = -6, v = 7,
drift = f_x, diff = g_x, Methods = "Euler")
11
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

Figure 5: Simulation 100-samples of the random variable X
t
at time v
t
= 7 by a simulated
100 trajectories of dX
t
= 2(1 X
t
)dt + dW
t
, used Euler algorithm, with X
0
= 6.
R> X
[1] 0.42457787 0.89022987 0.53798985 0.85601970 0.87574259
[6] 0.13023321 0.78466182 1.67646649 1.86814912 0.91651605
[11] 0.40458587 1.25308171 0.50399329 1.02546173 1.02647338
[16] 0.26630316 1.15195695 0.66666422 0.85174517 0.76239199
[21] 0.44828044 0.27588758 1.01659444 0.99617363 1.00829577
[26] 0.86676310 0.64266320 0.26960628 0.98420215 0.67478779
[31] 0.87145246 1.27619465 0.77522554 0.94704303 1.42153343
[36] 0.98755623 0.77789867 0.08774791 1.16267498 0.84184574
[41] 1.03631156 0.57812454 0.88114625 1.11535893 1.40388370
[46] 1.16252625 0.87735322 1.55669661 1.30152170 1.02433169
[51] 1.23616081 1.51367164 1.61080710 0.99608594 0.59199640
[56] 1.35987961 1.81565309 0.93158510 1.24933297 0.94759695
[61] 0.90637034 0.28367297 1.06232243 0.39463121 0.95167687
[66] 0.68991567 0.52513906 0.82081028 1.11318080 1.25055621
[71] 0.76350645 1.29918398 0.32224426 0.75747775 1.05616919
[76] 1.13505537 1.37253188 0.84595736 1.23985729 1.06338587
[81] 1.20910711 1.54841177 0.94479993 1.52138324 0.93499736
[86] 1.50113064 0.43927942 -0.09752953 0.99846875 1.69256126
12
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

[91] -0.24466704 1.31776649 0.77478081 0.65812555 0.98086811
[96] 1.34704883 1.24507068 -0.27647792 0.92965613 0.67631369
R> summary(X)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.2765 0.6865 0.9459 0.9235 1.2160 1.8680
Estimate parameters of the normal distribution by the method of maximum likelihood.
The function Ajdnorm needs as input the random variable X, and Initial values for opti-
mizer (mean, sd) with the condence level required. two methods are used of estimation
for the stationary distribution, the histograms used the function hist_general, and the
kernel methods by used Kern_general.
hist_general(Data, Breaks, Law = c("exp", "GAmma", "chisq", "Beta",
"fisher", "student", "weibull", "Normlog", "Norm"))
Kern_general(Data, bw, k, Law = c("exp", "GAmma", "chisq", "Beta",
"fisher", "student", "weibull", "Normlog", "Norm"))
R> Ajdnorm(X, starts = list(mean = 1, sd = 1), leve = 0.95)
Profiling...
$summary
Maximum likelihood estimation
Call:
mle(minuslogl = lik, start = starts)
Coefficients:
Estimate Std. Error
mean 0.9285440 0.04218993
sd 0.4197845 0.02983209
-2 log L: 109.0828
$coef
mean sd
0.9285440 0.4197845
$AIC
[1] 113.0828
$vcov
mean sd
13
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

mean 1.779990e-03 -1.799151e-10
sd -1.799151e-10 8.899539e-04
$confint
2.5 % 97.5 %
mean 0.8450425 1.0120456
sd 0.3674690 0.4858145
R> hist_general(Data = X, Breaks = 'Sturges', Law = "Norm")
R> Kern_general(Data = X, bw = "Bcv", k="gaussian", Law = "Norm")
Figure 6: Estimation stationary distribution used histograms and kernel methods.
3 Attractive models
The problem of dispersion is a very complex phenomenon is many problems dealing with
environment, biology, physics, chemistry, etc . . . , the dynamical behavior of such phe-
nomenon is a random process, often hard to modeling mathematically. This problem,
have been proposed by many authors [17, 19, 18, 3, 5]. For many dispersal problems, the
diusion processes are used to modeling the behavior of the dispersal phenomenon. For
example, two important models can be indicated, the rst is proposed by [5], and describes
the dispersion of a pollutant on an area of shallow water (attractive model for one diusion
process), the second model is the two diusions processes in attraction proposed by [6]
(attractive model for two diusions processes), which describes the dynamical behavior of
a two insects, one attracts the other.
14
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

3.1 Attractive model for one diusion process
Consider a shallow water area with depth L(x, y, z, t), horizontal U
w
(x, y, z, t) and V
w
(x, y, z, t),
S
w
(x, y, z, t) the velocities of the water in respectively the x, y and z directions, and
U
a
(x, y, z), V
a
(x, y, z), S
a
(x, y, z) the velocities of a particle caused by an attractive mech-
anism. Let (X
t
, Y
t
, Z
t
) be the position of a particle injected in the water at time t = t
0
at
the point (x
0
, y
0
, z
0
).
For a single particle, we propose the following dispersion models family [5]:
_

_
dX
t
=
_
U
a
+ U
w
+
L
x
L
D +
D
x
_
dt +

2DdW
1
t
dY
t
=
_
V
a
+ V
w
+
L
y
L
D +
D
y
_
dt +

2DdW
2
t
, t [0, T]
dZ
t
=
_
S
a
+ S
w
+
L
z
L
D +
D
z
_
dt +

2DdW
3
t
(10)
with:
U
a
=
Kx
_
_
x
2
+ y
2
+ z
2
_
s+1
, V
a
=
Ky
_
_
x
2
+ y
2
+ z
2
_
s+1
, S
a
=
Kz
_
_
x
2
+ y
2
+ z
2
_
s+1
.
and s 1, K > 0, W
1
t
and W
2
t
,W
3
t
are Brownian motions.
U
w
(x, y, z, t) and V
w
(x, y, z, t), S
w
(x, y, z, t) are neglected and the dispersion coecient
D(x, y, z) is supposed constant and equal to
1
2

2
, ( > 0).
Using Itos transform for (10), it is shown that the radial process R
t
= (X
t
, Y
t
, Z
t
) is
a Markovian diusion, solution of the stochastic dierential equation, given by:
dR
t
=
_

2
R
s1
t
2
K
R
s
t
_
dt + d

W
t
, t [0, T] (11)
where 2K >
2
and . is the Euclidean norm and

W
t
is a determined Brownian motion.
We take interest in the random variable
(s)
c
rst passage time (See [7, 8]) of the particle
through a sphere of radius c, centered at the origin of R
3
-space. The random variable
(s)
c
is dened by :

(s)
c
= inf(t 0|R
t
c and R
0
= r) (12)
This variable plays an important role in the prediction of the rate of particles reaching the
attractive center. We denote by M
s,
the family of models dened by (10) or (11).
3.1.1 Code example
For example the simulation of the model M
(s=1,)
can be made by discretization of equa-
tions (10) or (11) with s = 1, a simulate and estimate the density function of the random
variable
(1)
c
= 1/
(1)
c
(12).
Simulation the models (10) by used the function RadialP3D_1 with s = 1.
15
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

R> RadialP3D_1(N = 5000, t0 = 0, Dt = 0.001, X0 = 1, Y0 = 0.5, Z0 = 0.5,
v = 0.2, K = 2, Sigma = 0.2)
Figure 7: Simulation 3-dimensional attractive model M
(s=1,=0.2)
for one diusion process.
Simulation 100-sample for the rst passage time FPT of the model M
(s=1,)
used the
function tho_M1, and estimate the density function of the
(1)
c
by the Gamma law.
R> tho_M1(N = 1000, M = 100, t0 = 0, T = 1, R0 = 1, v = 0.2, K = 4, sigma =0.9,
Methods = "Euler")
R> FPT
[1] 0.071 0.124 0.046 0.117 0.165 0.126 0.178 0.059 0.218 0.118
[11] 0.305 0.064 0.135 0.188 0.271 0.161 0.099 0.123 0.086 0.137
[21] 0.135 0.101 0.135 0.265 0.204 0.084 0.139 0.186 0.133 0.099
[31] 0.340 0.099 0.133 0.073 0.160 0.217 0.091 0.181 0.139 0.097
[41] 0.118 0.235 0.133 0.371 0.096 0.301 0.106 0.206 0.137 0.233
[51] 0.082 0.156 0.135 0.154 0.060 0.113 0.077 0.132 0.230 0.135
[61] 0.212 0.169 0.238 0.236 0.211 0.098 0.176 0.117 0.205 0.176
[71] 0.220 0.116 0.167 0.073 0.220 0.160 0.099 0.169 0.102 0.197
[81] 0.128 0.172 0.080 0.196 0.120 0.219 0.161 0.283 0.254 0.113
[91] 0.392 0.050 0.138 0.111 0.074 0.216 0.202 0.098 0.091 0.056
R> Ajdgamma(X = 1/FPT, starts = list(shape = 1, rate = 1), leve = 0.95)
Profiling...
$summary
Maximum likelihood estimation
16
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

Call:
mle(minuslogl = lik, start = starts)
Coefficients:
Estimate Std. Error
shape 5.0320084 0.69279984
rate 0.6426731 0.09304918
-2 log L: 514.6746
$coef
shape rate
5.0320084 0.6426731
$AIC
[1] 518.6746
$vcov
shape rate
shape 0.4799716 0.06130030
rate 0.0613003 0.00865815
$confint
2.5 % 97.5 %
shape 3.7969651 6.5185981
rate 0.4768285 0.8423622
R> hist_general(Data = 1/FPT, Breaks = 'Sturges', Law = "GAmma")
R> Kern_general(Data = 1/FPT , bw ='Ucv', k ="gaussian", Law = "GAmma")
Figure 8: Estimation the density function of
(1)
c=0.2
by the Gamma law, used histograms
and kernel methods.
17
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

3.2 Attractive model for two diusions processes
In the following paragraph, we will propose a model of two diusions in attraction M

(.)
(V
(1)
t
)
and M

0
(V
(2)
t
), can modeling the behavior of two insects, one attracts the other.
Considers V
(1)
t
= (X
t,1
, X
t,2
, X
t,3
) and V
(2)
t
= (Y
t,1
, Y
t,2
, Y
t,3
) two random processes of dif-
fusion, which one supposes respectively representing the positions of a male insect and
an insect female, and the male is attracted by the female. The behavior of the female is
supposed to be a process of Brownian motion, dened by the following equation
dV
(2)
t
= I
33
dW
t
(13)
Whereas the behavior of the male is supposed to be a process of diusion, whose drift
is directed, at every moment t, towards the position of the female, and who is given by [6]
dV
(1)
t
= dV
(2)
t
+
m+1
(D
t
)D
t
dt + I
33
d

W
t
, (14)
where W
t
and

W
t
are two Brownian motion independent, and
_
D
t
= V
(1)
t
V
(2)
t

m
(d) =
K
d
m
(15)
K and m are positive constants.
The model suggested is following form
dV
(1)
t

_

_
dX
t,1
= dY
t,1

K(X
t,1
Y
t,1
)

(X
t,1
Y
t,1
)
2
+(X
t,2
Y
t,2
)
2
+(X
t,3
Y
t,3
)
2

m+1
dt + d

W
1
t
dX
t,2
= dY
t,2

K(X
t,2
Y
t,2
)

(X
t,1
Y
t,1
)
2
+(X
t,2
Y
t,2
)
2
+(X
t,3
Y
t,3
)
2

m+1
dt + d

W
2
t
dX
t,3
= dY
t,3

K(X
t,3
Y
t,3
)

(X
t,1
Y
t,1
)
2
+(X
t,2
Y
t,2
)
2
+(X
t,3
Y
t,3
)
2

m+1
dt + d

W
3
t
(16)
dV
(2)
t

_

_
dY
t,1
= dW
1
t
dY
t,2
= dW
2
t
dY
t,3
= dW
3
t
(17)
Using It os transform for (16) and (17), it is shown that the process X
t
= (V
(1)
t
, V
(2)
t
)
is a Markovian diusion, solution of the stochastic dierential equation, given by:
dX
t
=

2
X
m1
t
K
X
m
t
dt + dW
t
, t [0, T] (18)
where K >
2
and . is the Euclidean norm and W
t
is a determined Brownian motion.
18
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

This model makes it possible to carry out a dynamic simulation of the real phenomenon.
Using these simulations, one can also estimate the density of probability of the moment of
the rst meeting (V
(1)
t
, V
(2)
t
) between the two insects, dened by
(V
(1)
t
, V
(2)
t
) = lim
c0
inf{t 0|D
t
c}
3.2.1 Code example
For example the simulation two dimensional of the phenomenon M

(.)
(V
(1)
t
) M

0
(V
(2)
t
),
can be made by discretization of equations (16) and (17), by used the function TowDiffAtra2D
R> TowDiffAtra2D(N = 5000, t0 = 0, Dt = 0.001, T = 1, X1_0 = 2,
X2_0 = 2, Y1_0 = -0.5 , Y2_0 = -1, v = 0.05,
K = 3, m = 0.1, Sigma = 0.3)
Figure 9: Illustration of simulation two dimensional of a trajectory of the interaction
between two insects.
Simulation three dimensional of the phenomenon by used the function TowDiffAtra3D,
and simulation 100-sample for the moment of the rst meeting (V
(1)
t
, V
(2)
t
) between the
two insects used the function tho_02diff.
R> TowDiffAtra3D(N = 5000, t0 = 0, Dt = 0.001, T = 1, X1_0 = 1,
X2_0 = 0.5, X3_0 = 0, Y1_0 = -0.5, Y2_0 = 0.5,
Y3_0 = -1, v = 0.05, K = 3, m = 0.1, Sigma = 0.15)
19
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

Figure 10: Illustration of simulation three dimensional of a trajectory of the interaction
between two insects.
R> tho_02diff(N = 1000, M = 100, t0 = 0, Dt = 0.001, T = 1,
X1_0 = 1, X2_0 = 1, Y1_0 = 0.5, Y2_0 = 0.5,
v = 0.05, K = 4, m = 0.2, Sigma = 0.2)
R> FPT
[1] 0.143 0.109 0.123 0.106 0.133 0.123 0.189 0.120 0.270 0.140 0.198
[12] 0.144 0.107 0.092 0.165 0.178 0.124 0.146 0.113 0.136 0.158 0.169
[23] 0.108 0.160 0.143 0.201 0.122 0.091 0.154 0.096 0.189 0.198 0.147
[34] 0.147 0.128 0.105 0.192 0.106 0.139 0.174 0.134 0.105 0.090 0.165
[45] 0.284 0.098 0.136 0.092 0.093 0.077 0.149 0.171 0.125 0.151 0.122
[66] 0.221 0.199 0.154 0.140 0.145 0.217 0.106 0.097 0.121 0.131 0.153
[77] 0.152 0.168 0.127 0.100 0.120 0.130 0.181 0.175 0.166 0.102 0.136
[88] 0.127 0.136 0.137 0.081 0.112 0.204 0.113 0.109 0.185 0.130 0.186
[99] 0.193 0.140 0.106 0.191 0.131 0.203 0.123 0.128 0.132 0.158 0.126
R> Ajdgamma(X = 1/FPT, starts = list(shape = 1, rate = 1), leve = 0.95)
Profiling...
$summary
Maximum likelihood estimation
Call:
mle(minuslogl = lik, start = starts)
Coefficients:
20
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

Estimate Std. Error
shape 15.525968 2.1834585
rate 2.077233 0.2968927
-2 log L: 403.2067
$coef
shape rate
15.525968 2.077233
$AIC
[1] 407.2067
$vcov
shape rate
shape 4.7674910 0.63784671
rate 0.6378467 0.08814528
$confint
2.5 % 97.5 %
shape 11.657589 20.252660
rate 1.551179 2.719831
R> hist_general(Data = 1/FPT, Breaks = 'Sturges', Law = "GAmma")
R> Kern_general(Data = 1/FPT , bw ='Ucv', k ="gaussian", Law = "GAmma")
Figure 11: Estimation the density of probability of the moment of the rst meeting between
the two insects, used histograms and kernel methods.
21
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

4 Graphical User Interface for Sim.DiProc
Unlike S-PLUS, R
2
does not incorporate a statistical graphical user interface (GUI), but it
does include tools for building GUIs.
3
Based on the tcltk package (which furnishes an in-
terface to the Tcl/Tk GUI toolkit), the Sim.DiProcGUI package [10] provides a graphical
user interface for some functions in the Sim.DiProc package, the design objectives of the
Sim.DiProcGUI were as follows: to support Sim.DiProc, through an easy-to-use, exten-
sible, crossplatform GUI, keep things relatively simple, and to render visible, in a reusable
form. The Sim.DiProcGUI package uses a simple and familiar menu/dialog-box inter-
face. Top-level menus include File, Edit, Brownian Motion, Stochastic Integral, Stochastic
Models, Parametric Estimation, Numerical Solution of SDE, Statistical Analysis, and Help.
Each dialog box includes a Help button, which leads to a relevant help page, the R-
Sim.DiProc console also provides the ability to edit, enter, and re-execute commands.
Data sets in this GUI are simply R data frames, and can be read from attached packages
or imported from les, although several data frames may reside in memory.
Once R is running, simply loading the Sim.DiProcGUI package by typing the com-
mand library(Sim.DiProcGUI) into the R Console starts the GUI. After loading the
package, the GUI window should appear more or less as in the following gure 12
Figure 12: Graphical User Interface for Sim.DiProc package at start-up.
2
R is a programming language and software environment for statistical computing and graphics, are
available for download from CRAN at URL: http://www.r-project.org/.
3
The Sim.DiProcGUI package, described in this paper, is based on the tcltk package [12, 13], which
provides an interface to Tcl/Tk [26].
22
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

R> library(Sim.DiffProcGUI)
This and other screen images were created under Windows seven, if you use another
version of Windows (or, of course, another computing platform), then the appearance of
the screen may dier.
4
The Sim.DiProcGUI package and R Console windows oat freely on the desktop. You
will normally use the menus and dialog boxes of the Sim.DiProcGUI to read, manipulate,
and simulated or analyze data.
5 Conclusion
This paper introduces new package Sim.DiProc for a simulation of diusion processes in
R language, and graphical user interface (GUI) for this package, the actual development
of computing tools (software and hardware) has motivated us to analyze by simulation.
Many theoretical problems on the stochastic dierential equations have become the object
of practical research, as statistical analysis and the simulation the solution of SDE, enabled
many searchers in dierent domains to use these equations to modeling and to analyse
practical problems. The dispersion problem that we have treated in this paper is a good
example which shows the important use of the SDE in the practice, this problem is very
hard, hence the SDE approach seems to be a good approximation to treat such problem,
the diculty to obtain the exact solution of the SDE, the simulation gives information on
the density function of the random variable rst passage time
c
and enables to estimate a
density function, this estimation is based, on either, the histograms and the kernel methods.
These results can be applied, as a rst approximation, of the dispersal phenomenon, in
presence of an attractive source. The density function of
c
can be used to determine the
rate of the pollutant particles in a neighborhood of the attractive centre, it has been shown
by graphical and numerical simulations. The simulation studies implemented in R language
seems very preferment and ecient, because it is a statistical environment, which permits
to realize and to validate the simulations.
Acknowledgments
The work described in this paper was supported by the project entitled by Statistical
and Stochastic Modeling of the Complex Structures with code B00220080009, within the
framework National Commission of Evaluation and Futurology of the University Research
(CNEPRU), Ministry for Higher Education and Scientic Research, Algeria.
4
The examples in this document use the Windows version of R, R, however, is available on other
computing platforms as well (Macintosh computers and Unix/Linux systems), and the use of R and the
Sim.DiProcGUI package on these other systems is very similar to their use under Windows.
23
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

References
[1] E. Allen, Modeling with Ito Stochastic Dierential Equations, Springer-Verlag, New
York, 2007.
[2] K. Boukhetala, Estimation des param`etres de derives dune
`
Equation dierentielle
stochastique lineaire, controlee, masters thesis, University of Science and Technology
Houari Boumediene, BP 32 El-Alia, U.S.T.H.B, Algeria, 1984.
[3] , Simulation study of a dispersion about an attractive center, In proceeding of
11th Symposium Computational Statistics., (1994), pp. 128130. Edited by R. Dutter
and W. Grossman, Wien, Austria.
[4] , Identication and simulation of a communication system, Maghreb Mathemat-
ical Review, 2 (1995), pp. 5579.
[5] , Computer Methods and Water Resources, Modelling and Simulation of a Disper-
sion Pollutant with Attractive Centre, vol. 3, Boston, USA, computational mechanics
publications ed., 1996, pp. 245252.
[6] , Application des Processus de Diusion, Echantillonnage Optimal, PhD thesis,
University of Science and Technology Houari Boumediene, BP 32 El-Alia, U.S.T.H.B,
Algeria, 1998.
[7] , Estimation of the rst passage time distribution for a simulated diusion process,
Maghreb Mathematical Review, 7 (1998), pp. 125.
[8] , Kernel density of the exit time in a simulated diusion, The Annals of The
Engineer Maghrebian, 12 (1998), pp. 587589.
[9] K. Boukhetala and A. Guidoum, Sim.DiProc: Simulation of Diusion Pro-
cesses, 2011. R package version 2.0.
[10] , Sim.DiProcGUI: Graphical User Interface for Simulation of Diusion Pro-
cesses, 2011. R package version 2.0.
[11] F. Brodeau and A. L. Breton, Identication de param`etres pour un syst`eme excite
par des bruits gaussien et poissonien, The Annals of The Institute Henri Poincare, 1
(1979), pp. 123.
[12] P. Dalgaard, A primer on the R-tcl/tk package, R News, 1 (2001), pp. 2731.
[13] , Changes to the R-tcl/tk package, R News, 2 (2002), pp. 2571.
[14] H. Douglas and P. Peter, Stochastic Dierential Equations in Science and En-
gineering, World Scientic Publishing, 2006.
24
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-

[15] J. Franck, Mod`eles aleatoires et physique probabiliste, Springer-Verlag, New York,
2009.
[16] A. Greiner, W. Strittmatter, and J. Honerkamp, Numerical integration of
stochastic dierential equations, The Journal of Statistical Physics, 51 (1988).
[17] K. Hadeler, P. de Mottoni, and K. Schumacher, Dynamic models for animal
orientation, The Journal of Mathematical Biology, 10 (1980), pp. 307332.
[18] A. Heemink, Stochastic modelling of dispersion in shallow water, Stochastic Hydrol-
ogy and Hydraulics, 4 (1990), pp. 161174.
[19] S. Helland, Diusion models for the dispersal of insects near an attractive center,
The Journal of Mathematical Biology, 18 (1983), pp. 103122.
[20] P. Kloeden and E. Platen, A survey of numerical methods for stochastic dier-
ential equations, Stochastic Hydrology and Hydraulics, 3 (1989), pp. 155178.
[21] E. Peter and P. Eckhard, Numerical Solution of Stochastic Dierential Equa-
tions, Springer-Verlag, New York, 1995.
[22] R Development Core Team, R: A Language and Environment for Statistical
Computing, R Foundation for Statistical Computing, Vienna, Austria, 2011. ISBN
3-900051-07-0.
[23] T. Rolski, H. Schmidli, V. Schmidt, and J. Teugels, Stochastic Processes for
Insurance and Finance, John Wiley & Sons, 1998.
[24] Y. Saito and T. Mitsui, Simulation of stochastic dierential equations, The Annals
of the Institute of Statistical Mathematics, 3 (1993), pp. 419432.
[25] M. Stefano, Simulation and Inference for Stochastic Dierential Equations,
Springer-Verlag, New York, 2008.
[26] B. Welch, Practical Programming in Tcl and Tk, Prentice Hall PTR Upper Saddle
River, NJ, USA, 3 ed., 2000.
25
h
a
l
-
0
0
6
2
9
8
4
1
,

v
e
r
s
i
o
n

1

-