sde.sim {sde}R Documentation

Simulation of Stochastic Differential Equation

Description

Generic interface to different methods of simulations of solutions to stochastic differential equations.

Usage

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)

Arguments

t0 time origin
T horizon of simulation
X0 initial value of the process
N number of simulation steps
delta time-step of the simulation
drift drift coeffcient: a expression of two variables t and x
sigma diffusion coeffcient: a expression of two variables t and x
drift.x partial derivative of drift coeffcient wrt to x: a function of two variables t and x
sigma.x partial derivative of diffusio coefficient wrt to x: a function of two variables t and x
drift.xx second partial derivative of drift coefficient wrt to x: a function of two variables t and x
sigma.xx second partial derivative of diffusion coefficient wrt to x: a function of two variables t and x
drift.t partial derivative of drift coefficient wrt to t: a function of two variables t and x
method method of simulation, see details
alpha weight alpha of the predictor-corrector scheme
eta weight eta of the predictor-corrector scheme
pred.corr boolean: wheter to apply the predictor-correct adjustment. See details.
rcdist a function which is a random number generator from the conditional distribution of the process. See details.
theta vector of parameters for cdist. See details.
model model from which to simulate. See details.
k1 lower bound for psi(x) = 0.5*drift(x)^2 + 0.5*drift.x(x). See details.
k2 upper bound for psi(x) = 0.5*drift(x)^2 + 0.5*drift.x(x). See details.
phi the function psi(x) - k1
max.psi upper value of the support of psi to search for its maximum
rh the rejection function. Seed details.
A A(x) is the integral of the drift between 0 and x

Details

The function returns a ts object of length N+1, i.e. X0 and the new N simulated values. If delta is not specified, then delta = (T-t0)/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.

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 the constant function 1

The method of simulation can be one among: euler, milstein, milstein2, KPS, cdist, ozaki and shoji. No assumption on the coefficients or on cdist is checked: the user is responbile for using the right method for the process object of simulation.

The model is one among: CIR: Cox-Ingersoll-Ross, VAS: Vasicek, OU Ornstein-Uhlenbeck, BS: Black and Scholes. No assumption on the coefficients theta is checked: the user is responbile for using the right ones.

If 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 exact algorithm method EA: if missing k1 and k2 as well as A, rh and phi are calculated numerically by the function.

Value

x returns and invisible ts object

Note

This package is a companion to the book `Simulation and Inference for Stochastic Differential Equation, Springer, NY.

Author(s)

Stefano Maria Iacus

References

See Chapter 2 of the book

Examples

# Ornstein-Uhlenbeck process
set.seed(123)
d <- expression(-5 * x)
s <- expression(3.5) 
sde.sim(X0=10,drift=d, sigma=s) -> X
plot(X,main="Ornstein-Uhlenbeck")

# Cox-Ingersoll-Ross process
# dXt = (6-3*Xt)*dt + 2*sqrt(Xt)*dWt
set.seed(123)
d <- expression( 6-3*x ) 
s <- expression( 2*sqrt(x) ) 
sde.sim(X0=10,drift=d, sigma=s) -> X
plot(X,main="Cox-Ingersoll-Ross")

# Cox-Ingersoll-Ross 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="Cox-Ingersoll-Ross")

set.seed(123)
sde.sim(X0=10, theta=c(6, 3, 2), model="CIR") -> X
plot(X, main="Cox-Ingersoll-Ross")

# Exact simulation
set.seed(123)
d <- expression(sin(x))
d.x <- expression(cos(x)) 
A <- function(x) 1-cos(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")

[Package sde version 1.8 Index]