Generic interface to different methods of simulation of solutions to stochastic differential equations.
1 2 3 4 5 6 7  sde.sim(t0 = 0, T = 1, X0 = 1, N = 100, delta, drift, sigma,
drift.x, sigma.x, drift.xx, sigma.xx, drift.t,
method = c("euler", "milstein", "KPS", "milstein2",
"cdist","ozaki","shoji","EA"),
alpha = 0.5, eta = 0.5, pred.corr = T, rcdist = NULL,
theta = NULL, model = c("CIR", "VAS", "OU", "BS"),
k1, k2, phi, max.psi = 1000, rh, A, M=1)

t0 
time origin. 
T 
horizon of simulation. 
X0 
initial value of the process. 
N 
number of simulation steps. 
M 
number of trajectories. 
delta 
time step of the simulation. 
drift 
drift coefficient: an expression of two variables 
sigma 
diffusion coefficient: an expression of two variables 
drift.x 
partial derivative of the drift coefficient w.r.t. 
sigma.x 
partial derivative of the diffusion coefficient w.r.t. 
drift.xx 
second partial derivative of the drift coefficient w.r.t. 
sigma.xx 
second partial derivative of the diffusion coefficient w.r.t. 
drift.t 
partial derivative of the drift coefficient w.r.t. 
method 
method of simulation; see details. 
alpha 
weight 
eta 
weight 
pred.corr 
boolean: whether to apply the predictorcorrect adjustment; see details. 
rcdist 
a function that is a random number generator from the conditional distribution of the process; see details. 
theta 
vector of parameters for 
model 
model from which to simulate; see details. 
k1 
lower bound for 
k2 
upper bound for 
phi 
the function 
max.psi 
upper value of the support of 
rh 
the rejection function; see details. 
A 

The function returns a ts
object of length N+1
; i.e., X0
and
the new N
simulated values if M=1
.
For M>1
, an mts
(multidimensional ts
object) is returned, which
means that M
independent trajectories are simulated.
If the initial value X0
is not of the length M
, the values are recycled
in order to have an initial vector of the correct length.
If delta
is not specified, then delta = (Tt0)/N
.
If delta
is specified, then N
values of the solution of the sde are generated and
the time horizon T
is adjusted to be N * delta
.
The function psi
is psi(x) = 0.5*drift(x)^2 + 0.5*drift.x(x)
.
If any of drift.x
, drift.xx
, drift.t
,
sigma.x
, and sigma.xx
are not specified,
then numerical derivation is attempted when needed.
If sigma
is not specified, it is assumed to be the constant function 1
.
The method
of simulation can be one among: euler
, KPS
, milstein
,
milstein2
, cdist
, EA
, ozaki
, and shoji
.
No assumption on the coefficients or on cdist
is checked: the user is
responsible for using the right method for the process object of simulation.
The model
is one among: CIR
: CoxIngersollRoss, VAS
: Vasicek,
OU
OrnsteinUhlenbeck, BS
: Black and Scholes.
No assumption on the coefficient theta
is checked: the user is responsible
for using the right ones.
If the method
is cdist
, then the process is simulated according to its
known conditional distribution. The random generator rcdist
must be a
function of n
, the number of random numbers; dt
, the time lag;
x
, the value of the process at time t
 dt
; and the
vector of parameters theta
.
For the exact algorithm method EA
: if missing k1
and k2
as well
as A
, rh
and phi
are calculated numerically by the function.
x 
returns an invisible 
Stefano Maria Iacus
See Chapter 2 of the text.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39  # OrnsteinUhlenbeck process
set.seed(123)
d < expression(5 * x)
s < expression(3.5)
sde.sim(X0=10,drift=d, sigma=s) > X
plot(X,main="OrnsteinUhlenbeck")
# Multiple trajectories of the OU process
set.seed(123)
sde.sim(X0=10,drift=d, sigma=s, M=3) > X
plot(X,main="Multiple trajectories of OU")
# CoxIngersollRoss process
# dXt = (63*Xt)*dt + 2*sqrt(Xt)*dWt
set.seed(123)
d < expression( 63*x )
s < expression( 2*sqrt(x) )
sde.sim(X0=10,drift=d, sigma=s) > X
plot(X,main="CoxIngersollRoss")
# CoxIngersollRoss using the conditional distribution "rcCIR"
set.seed(123)
sde.sim(X0=10, theta=c(6, 3, 2), rcdist=rcCIR,
method="cdist") > X
plot(X, main="CoxIngersollRoss")
set.seed(123)
sde.sim(X0=10, theta=c(6, 3, 2), model="CIR") > X
plot(X, main="CoxIngersollRoss")
# Exact simulation
set.seed(123)
d < expression(sin(x))
d.x < expression(cos(x))
A < function(x) 1cos(x)
sde.sim(method="EA", delta=1/20, X0=0, N=500,
drift=d, drift.x = d.x, A=A) > X
plot(X, main="Periodic drift")

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.